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Abstract 

We  develop  divergence-conforming  B-spline  discretizations  for  the  numerical  solu¬ 
tion  of  the  steady  Navier-Stokes  equations.  These  discretizations  are  motivated  by 
the  recent  theory  of  isogeometric  discrete  differential  forms  and  may  be  interpreted 
as  smooth  generalizations  of  Raviart-Thomas  elements.  They  are  (at  least)  patch- 
wise  and  can  be  directly  utilized  in  the  Galerkin  solution  of  steady  Navier-Stokes 
flow  for  single-patch  conhgurations.  When  applied  to  incompressible  flows,  these  dis¬ 
cretizations  produce  pointwise  divergence- free  velocity  helds  and  hence  exactly  satisfy 
mass  conservation.  Consequently,  discrete  variational  formulations  employing  the  new 
discretization  scheme  are  automatically  momentum-conservative  and  energy-stable. 

In  the  presence  of  no-slip  boundary  conditions  and  multi-patch  geometries,  the  dis¬ 
continuous  Galerkin  framework  is  invoked  to  enforce  tangential  continuity  without 
upsetting  the  conservation  or  stability  properties  of  the  method  across  patch  bound¬ 
aries.  Furthermore,  as  no-slip  boundary  conditions  are  enforced  weakly,  the  method 
automatically  defaults  to  a  compatible  discretization  of  Euler  flow  in  the  limit  of 
vanishing  viscosity.  The  proposed  discretizations  are  extended  to  general  mapped  ge¬ 
ometries  using  divergence-preserving  transformations.  For  sufficiently  regular  single¬ 
patch  solutions  subject  to  a  smallness  condition,  we  prove  a  priori  error  estimates 
which  are  optimal  for  the  discrete  velocity  field  and  suboptimal,  by  one  order,  for  the 
discrete  pressure  held.  We  present  a  comprehensive  suite  of  numerical  experiments 
which  indicate  optimal  convergence  rates  for  both  the  discrete  velocity  and  pressure 
helds  for  general  conhgurations,  suggesting  that  our  a  priori  estimates  may  be  conser¬ 
vative.  These  numerical  experiments  also  suggest  our  discretization  methodology  is 
robust  with  respect  to  Reynolds  number  and  more  accurate  than  classical  numerical 
methods  for  the  steady  Navier-Stokes  equations. 

Keywords:  Steady  Navier-Stokes  equations,  B-splines,  Isogeometric  analysis.  Divergence- 
conforming  discretizations 
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1  Introduction 


Steady  Navier-Stokes  flow  is  an  important  simplification  of  fnlly  unsteady  Navier- 
Stokes  flow.  Many  low  speed,  laminar  fluid  flows  may  be  accurately  described  by 
the  steady  Navier-Stokes  equations.  Additionally,  one  arrives  at  a  steady  Navier- 
Stokes  problem  by  conducting  a  Reynolds  time-averaging  of  statistically  stationary 
Navier-Stokes  flows  (see  Chapters  3  and  4  of  [31]).  Despite  its  simple  appearance, 
the  steady  Navier-Stokes  problem  has  presented  considerable  difficulty  in  its  numer¬ 
ical  approximation.  It  is  subject  to  the  Babuska-Brezzi  inf-sup  condition,  and  when 
the  convection  operator  is  expressed  in  conservation  form  and  the  incompressibility 
constraint  is  not  met  exactly,  terms  corresponding  to  convection  can  actually  accrete 
energy.  This  ultimately  leads  to  unstable  numerical  formulations,  and  alternative  rep¬ 
resentations  of  the  convection  operator  have  been  devised  to  bypass  this  instability. 
The  most  popular  of  these  in  the  finite  element  community  is  the  skew-symmetric 
representation  [40].  Discretizations  of  the  convection  term  using  the  skew-symmetric 
representation  neither  produce  nor  dissipate  energy  and  hence  lead  to  stable  numer¬ 
ical  methods.  Unfortunately,  these  discretizations  do  not  inherit  the  conservation 
structure  of  the  Navier-Stokes  equations.  Alternatively,  provably  stable,  convergent, 
and  locally-conservative  discontinuous  Galerkin  discretizations  have  been  devised  for 
the  steady  Navier-Stokes  equations  in  [11,  12],  but  these  discretizations  are  encum¬ 
bered  with  a  proliferation  of  degrees  of  freedom  and  are  thus  largely  limited  to  two 
spatial  dimensions.  Hybrid  technologies  have  recently  been  proposed  with  the  aim  of 
extending  the  applicability  of  discontinuous  Galerkin  methods  to  larger  problem  sizes 
[28,  29]. 

Another  discretization  procedure  for  the  steady  Navier-Stokes  equations  arises 
through  the  intelligent  choice  of  weighting  function  in  a  Petrov-Galerkin  method.  A 
popular  method  of  choice  is  the  use  of  an  advective  formulation  in  conjunction  with 
the  Streamline-Upwind  Petrov-Galerkin  (SUPG)  method  [7]  to  handle  convective  in¬ 
stabilities  and  the  Pressure-Stabilizing  Petrov-Galerkin  (PSPG)  method  [25,  37]  to 
handle  pressure  instabilities.  Unfortunately,  the  theoretical  analysis  of  this  method 
in  the  steady  regime  has  been  entirely  restricted  to  linearized  Oseen  problems  where 
the  convection  velocity  is  assumed  fixed  and  divergence-free.  These  linearized  model 
problems  are  ultimately  insufficient  as  the  discrete  convection  velocity  is  not,  in  gen¬ 
eral,  divergence-free.  To  further  control  the  divergence  term,  so-called  grad-div  stabi¬ 
lization  techniques  [33]  have  been  proposed  which  add  artificial  dilatational  stresses 
to  the  underlying  variational  formulation.  Using  a  combination  of  an  advective  for¬ 
mulation,  SUPG,  PSPG,  and  grad-div  stabilization,  provably  convergent  numerical 
methods  have  been  devised  for  the  steady  Navier-Stokes  equations  (see  Ghapter  3 
of  [33]),  and  such  methods  can  be  made  to  globally  conserve  momentum  through  a 
residual-based  modification  [26].  Still,  a  provably  convergent  H^-Galerkin  finite  ele¬ 
ment  discretization  of  the  three-dimensional  steady  Navier-Stokes  equations  written 
in  conservation  form  has  proved  elusive. 
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In  this  paper,  we  present  divergence-conforming  B-spline  discretizations  for  the 
steady  Navier-Stokes  problem.  These  discretizations  are  motivated  by  the  theory  of 
isogeometric  discrete  differential  forms  [9,  10]  and  extend  the  Darcy-Stokes-Brinkman 
discretizations  presented  in  [19]  to  nonlinear  Navier-Stokes  flows.  As  our  discretiza¬ 
tions  return  pointwise  divergence-free  velocity  helds,  we  can  utilize  a  variational  for¬ 
mulation  written  in  conservation  form  without  being  susceptible  to  instability.  We 
impose  no-penetration  boundary  conditions  strongly  and  no-slip  boundary  conditions 
weakly  using  Nitsche’s  method.  This  allows  our  discretization  procedure  to  naturally 
default  to  a  conforming  approximation  of  Euler  flow  in  the  limit  of  vanishing  vis¬ 
cosity.  This  also  allows  our  method  to  capture  boundary  layers  without  resorting 
to  stretched  meshes  [3,  4].  We  prove  stability  and  error  estimates  for  single-patch 
discretizations  under  a  smallness  condition.  Our  error  estimates  are  optimal  for  the 
discrete  velocity  held  and  suboptimal,  by  one  order,  for  the  discrete  pressure  held 
provided  that  that  the  exact  solution  is  sufficiently  regular.  All  of  our  estimates’  de¬ 
pendencies  on  the  viscosity  and  the  penalty  parameter  of  Nitsche’s  method  are  made 
explicit  in  our  analysis.  We  utilize  the  methods  of  exact  and  manufactured  solutions 
to  verify  our  error  estimates  and  hnd  our  discrete  pressure  helds  converge  at  optimal 
order  in  contrast  with  our  theoretical  estimates.  We  further  test  the  ehectiveness  of 
our  method  by  considering  the  application  of  our  discretization  to  the  analysis  of  two 
benchmark  problems:  lid-driven  cavity  how  and  conhned  jet  impingement. 

An  outline  of  this  paper  is  as  follows.  In  the  following  section,  we  present  some 
basic  notation.  In  Section  3,  we  recall  the  steady  Navier-Stokes  problem  subject  to 
homogeneous  Dirichlet  boundary  conditions.  In  Section  4,  we  briehy  review  B-splines, 
the  basic  building  blocks  of  our  new  discretization  technique,  and  in  Section  5,  we 
dehne  the  B-spline  spaces  which  we  will  utilize  to  discretize  velocity  and  pressure 
helds.  In  Section  6,  we  present  our  discrete  variational  formulation  for  the  steady 
Navier-Stokes  problem  and  prove  continuity,  stability,  and  a  priori  error  estimates 
for  the  single-patch  setting.  In  Section  7,  we  discuss  the  extension  of  our  methodology 
to  the  multi-patch  setting.  In  Sections  8  and  9,  we  present  numerical  results,  and  in 
Section  10,  we  draw  conclusions.  Before  proceeding,  note  that  one  might  say  there  is 
a  fundamental  issue  concerning  the  fact  that  our  analysis  only  covers  hows  subject  to 
“small  data” .  However,  well-posedness  of  the  continuous  problem  is  subject  to  a  sim¬ 
ilar  constraint,  and  we  believe  the  small  data  assumption  is  natural  as  medium-  and 
large-Reynolds  number  hows  are  inherently  unsteady  in  both  laminar  and  turbulent 
regimes. 


2  Notation 

We  begin  this  paper  with  some  basic  notation.  For  d  a  positive  integer  representing 
dimension,  let  H  C  denote  an  arbitrary  bounded  Lipschitz  domain  with  boundary 
dD.  As  usual,  let  L^{D)  denote  the  space  of  square  integrable  functions  on  D  and 
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define  L^(Zi))  =  {L‘^{D)Y.  We  will  also  utilize  the  more  general  Lebesgue  spaces 
L^{D)  where  1  <  p  <  cx)  and  their  vectorial  counterparts  Ij^{D).  Let  H^{D)  denote 
the  space  of  functions  in  L‘^{D)  whose  k^^-order  derivatives  belong  to  L‘^{D)  and 
dehne  =  (yH^{D)Y ■  We  identify  with  H^{D)  the  standard  Sobolev  norm 

ll'^llrf''(D)  —  I  ^  II-D"v||^2(£,)  I 

yict|<fc  J 

where  cx  =  (ai,  0:2, . . . ,  (Xd)  is  a  multi-index  of  non-negative  integers,  |q:|  =  ai  +  a2  + 
. . .  +  ad,  and 

^l«l 

^  8x2'^ . . .  dx‘^‘^ ' 

We  denote  the  Sobolev  semi-norms  as  |  ■  \h>^{d),  and  we  adopt  the  convention  H^{D)  = 
L‘^(D).  Throughout,  Sobolev  spaces  of  fractional  order  are  dehned  using  function 
space  interpolation  (see,  e.g..  Chapter  1  of  [38]).  We  dehne  Hq{D)  C  H^{D)  to  be 
the  subspace  of  functions  with  homogeneous  boundary  conditions  and  dehne  Hj(Zi)) 
to  be  the  vectorial  counterpart  of  Hq{D).  We  dehne  H®(div;Zi))  to  be  the  Sobolev 
space  of  all  functions  in  H®(r2)  whose  divergence  also  belongs  to  H^{D).  This  space 
is  equipped  with  the  norm 

||v||Hydiv;D)  =  (llvllnyD)  +  l|divv||^«(^))  . 

When  s  =  0,  we  omit  the  superscript.  We  also  dehne 

Ho(div;  D)  =  {v  G  H(div;  D)  :  v  ■  n  =  0  on  dD} 

where  n  denotes  the  outward  pointing  unit  normal.  Finally,  we  denote  Lq{D)  C  L‘^{D) 
as  the  space  of  square-integrable  functions  with  zero  average  on  D. 

3  The  Steady  Navier- Stokes  Problem 

In  this  section,  we  recall  the  steady  Navier-Stokes  problem  subject  to  homogeneous 
Dirichlet  boundary  conditions.  For  d  a  positive  integer,  let  hi  denote  a  Lipschitz 
bounded  open  set  of  Throughout  this  paper,  d  will  be  either  2  or  3.  The  problem 
of  interest  is  as  follows. 

Given  1/  G  M"*"  and  f  ;  hi  — )■  hnd  u  :  0  — )■  and  p  :  hi  — )■  M  such 

that 

V  ■  (u  ®  u)  —  V  ■  (2z/V'^u)  -|-  gradp  =  f  in  hi  (1) 

divu  =  0  in  n  (2) 

u  =  0  on  dQ.  (3) 
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Above,  u  denotes  the  flow  velocity  of  a  fluid  moving  through  the  domain  hi,  p  de¬ 
notes  the  pressure  acting  on  the  fluid  divided  by  the  fluid  density,  u  denotes  the 
kinematic  viscosity  of  the  fluid,  f  denotes  a  body  force  acting  on  the  fluid  divided  by 
the  density,  and  V^u  denotes  the  symmetrized  gradient  of  the  velocity  held  defined 
by  V^u  =  ^  ^Vu-|-  (Vu)^j  .  Note  that  the  pressure  is  only  determined  up  to  an 
arbitrary  constant. 

Assuming  that  f  G  L^(r2) ,  the  weak  form  for  the  steady  Navier-Stokes  problem  is 
written  as  follows: 

(  Find  u  G  HQ(r2)  and  p  G  Ll{Q)  such  that 


A;(u, 

v)  c(u,  u;  v)  -  b{p,  v)  4- 

-b{q,u) 

=  (f,v; 

(4) 

for  all  V  G  Hg 

(12)  and  q  G  Fq(12)  where 

k(w,  v)  = 

(2z/V"w,VW)(^2(f^))dxd, 

Vw,  V 

(5) 

6(g,v)  = 

(g,  divv)^2(^) , 

yq  e  1 

1^(12),  V 

(6) 

c(w,x;  v)  = 

-(w®X, 

Vw,  X, 

,vGH1 

,(n). 

(7) 

V 


Note  that  the  trilinear  form  c(-,  ■;  ■)  ;  HQ(r2)  x  — )■  M  makes  sense  due 

to  the  continuous  Sobolev  embedding 


(8) 


In  fact,  as  dil  is  Lipschitz,  we  have  the  stronger  embedding  H^(r2)  L‘^(r2). 

We  have  the  following  existence  and  uniqueness  theorem  for  flows  subject  to  small 
data  whose  proof  may  be  found  in  [22], 

Theorem  3.1.  Problem  (W)  has  a  unique  weak  solution  (u,p)  G  Hg(r2)  X  mn) 
provided  the  problem  data  satisfies  an  inequality  of  the  form 


CnC 


Poin  I 


f||L2(0)  <  1 


(9) 


where  Cq  is  a  constant  which  only  depends  on  Vt  and  Cpoin  is  the  positive  constant 
appearing  in  Poincare’s  inequality: 


ll’'^llHho)  ^  I v|h1(o) 5  Vv  G  11^(12)  fl  Ho(div;  hi) . 

Furthermore,  such  a  weak  solution  satisfies  the  inequality 

Cpoin  1 1  p  1 1 


I^IhFo)  ^ 


(10) 


(11) 
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Remark  3.1.  Note  that,  since  v  is  constant  and  divu  =  0, 

V  ■  (2z/V'^u)  =  uAu.  (12) 

This  inspires  a  different  variational  formulation  than  that  presented  here  which  is 
often  the  basis  for  numerical  discretization  (see,  for  example,  [12]).  However,  such  a 
formulation  cannot  easily  accommodate  traction  boundary  conditions.  We  have  found 
that  discrete  formulations  based  on  either  of  the  diffusion  operators  given  by  (12)  yield 
qualitatively  and  quantitatively  similar  results. 

Remark  3.2.  In  general,  we  may  replace  the  constant-viscosity  Newtonian  stress 
tensor  given  here,  T  =  2z/V®u,  with  more  suitable  choices  of  stress  tensor.  Our 
analysis  does  not  cover  this  general  setting. 

4  B-splines  and  Geometrical  Mappings 

In  this  section,  we  briefly  introdnce  B-splines,  the  primary  ingredient  in  our  discretiza¬ 
tion  technique  for  the  steady  Navier-Stokes  equations.  We  also  introduce  mappings 
which  will  allow  us  to  extend  our  discretization  technique  to  general  geometries  of 
engineering  interest.  For  an  overview  of  B-splines,  their  properties,  and  robust  algo¬ 
rithms  for  evaluating  their  values  and  derivatives,  see  de  Boor  [14]  and  Schumaker 
[35].  For  the  application  of  B-splines  to  finite  element  analysis,  see  Hollig  [23]  and 
Cottrell,  Hughes,  and  Bazilevs  [13]. 

4.1  Univariate  B-splines 

For  two  positive  integers  k  and  n,  representing  degree  and  dimensionality  respectively, 
let  us  introduce  the  ordered  knot  vector 

S:={0  =  6,6,---,W+i  =  1}  (13) 


where 


■Cl  —  '^2  ^  •  'Cn+fc+l- 


Given  S  and  k,  univariate  B-spline  basis  functions  are  constructed  recursively  starting 
with  piecewise  constants  (/c  =  0): 


1  if  6  <  e  <  6+1 

0  otherwise. 


For  A;  =  1,  2,  3, . . .,  they  are  defined  by 


BfiO 


^  6  Dfc— l/'fl  I  6+A:+1  ^ 

(  (B^\h,)+  -°i+i  (s)- 

Si+fc  Si  Si+fc+1  Si+1 


(14) 


(15) 
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When  =  0,  ,  ^  \  is  taken  to  be  zero,  and  similarly,  when  fj+fc+i  —  6+1  =  0, 

Si  +  fe  si 

b+fc+i-g  jg  to  be  zero.  B-spline  basis  functions  are  piecewise  polynomials  of 

degree  k,  form  a  partition  of  unity,  have  local  support,  and  are  non-negative.  We 
refer  to  linear  combinations  of  B-spline  basis  functions  as  B-splines  or  simply  splines. 

Let  us  now  introduce  the  vector  <^  =  {^1, . . . ,  (m}  of  knots  without  repetitions  and 
a  corresponding  vector  {ri, . . . ,  r^}  of  knot  multiplicities.  That  is,  r*  is  dehned  to  be 
the  multiplicity  of  the  knot  Q  in  S.  We  assume  that  <  k  +  1.  Let  us  further  assume 
throughout  that  ri  =  =  /c  +  1,  he,  that  S  is  an  open  knot  vector.  This  allows 

us  to  easily  prescribe  Dirichlet  boundary  conditions.  At  the  point  Q,  B-spline  basis 
functions  have  aj  :=  k  —  rj  continuous  derivatives.  We  dehne  the  regularity  vector  ck 
by  CK  :=  {tti, . . .  ,am}-  By  construction,  ai  =  am  =  —1-  In  what  follows,  we  utilize 
the  notation 

|q:|  =  minjo;*  :  2  <  i  <  m  —  1}  (16) 

and  CK  —  1  :=  {—1,  0:2  —  1, ... ,  Qim-i  —  1,  —1}  when  to  >  0  for  2  <  i  <  m  —  1. 

We  denote  the  space  of  B-splines  spanned  by  the  basis  functions  as 

Sy=spati{B''to.  (17) 

When  k  >  1  and  a*  >  0  for  2  <  i  <  m  —  1,  the  derivatives  of  functions  in  are 
splines  as  well.  In  fact,  we  have  the  stronger  relationship 

One  of  the  most  important  properties  of  univariate  B-splines  is  rehnement  and,  per¬ 
haps  more  importantly,  nestedness  of  rehnement.  Knot  insertion  and  degree  elevation 
algorithms  are  described  in  detail  in  Chapter  2  of  [13]. 

4.2  Multivariate  B-splines 

The  dehnition  of  multivariate  B-splines  follows  easily  through  a  tensor-jproduct  con¬ 
struction.  For  d  again  a  positive  integer,  let  us  consider  the  unit  cube  O  =  (0, 1)'^  C 
which  we  will  henceforth  refer  to  as  the  parametric  domain.  Mimicking  the  one¬ 
dimensional  case,  given  integers  ki  and  n;  for  /  =  l,...,(i,  let  us  introduce  open 
knot  vectors  S;  =  . . . ,  ^ni+ki+i,i}  and  the  associated  vectors  Ci  =  •  •  • ,  Cm,,;}, 

{ri^i, . . . ,  rmi,i},  and  cc;  =  {cn,;, . . . ,  am,,;}-  There  is  a  parametric  Cartesian  mesh  Aih 
associated  with  these  knot  vectors  partitioning  the  parametric  domain  into  rectangu¬ 
lar  parallelepipeds.  Visually, 

M/i  =  {Q  =  (Ci,,;,  Ci,+i,;) ,  1  <  b  <  —  1}  •  (19) 

For  each  element  Q  €  Ad/i  we  associate  a  parametric  mesh  size  hq  =  diam((5).  We 

also  dehne  a  shape  regularity  constant  A  which  satishes  the  inequality 

A-i  <  <  A,  VQ  e  M,„  (20) 

hq 
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where  hg^min  denotes  the  length  of  the  smallest  edge  of  Q-  A  sequence  of  parametric 
meshes  that  satisfy  the  above  inequality  for  an  identical  shape  regularity  constant  is 
said  to  be  locally  quasi-uniform. 

We  associate  with  each  knot  vector  (/  =  1, ...  ,d)  univariate  B-spline  basis 
functions  of  degree  ki  for  ii  =  1, ...  ,ni.  On  the  mesh  we  dehne  the  tensor- 
product  B-spline  basis  functions  as 

Bfed;  ■= -St'.  ®  ■  ■  ■  ®  *1  =  1.. ..,n„  ...  ij  =  l,...,ni.  (21) 


We  then  accordingly  dehne  the  tensor-product  B-spline  space  as 


_  qki,...,ka 


{Mh) 


span 


■■,kd 

■,id 


ni,...,nd 

il=l,...,id=i 


(22) 


Like  their  univariate  counterparts,  multivariate  B-spline  basis  functions  are  piece- 
wise  polynomial,  form  a  partition  of  unity,  have  local  support,  and  are  non-negative. 
Dehning  the  regularity  constant 


a  :=  min  min  {aj, ;}  (23) 

we  see  that  our  B-splines  are  O^-continuous  throughout  the  domain  O.  Rehnement  of 
multivariate  B-spline  bases  is  obtained  by  applying  knot  insertion  and  degree  elevation 
in  tensor-product  fashion.  In  the  remainder  of  the  text,  we  consider  a  family  of 
nested  meshes  {Aih}h<ho  associated  B-spline  spaces 

have  been  obtained  by  successive  applications  of  knot  rehnement.  Furthermore,  we 
assume  throughout  that  the  mesh  family  {■Mh}h<ho  locally  quasi-uniform. 

Note  that  each  element  Q  =  <S)i=i,...,d  {Ch,^  Cii+i,i)  has  the  equivalent  representation 
Q  =  ®z=i,...,d  ^ji+i,i)  for  some  index  j).  With  this  in  mind,  we  associate  with  each 
element  a  support  extension  Q,  dehned  as 

Q  :=  {^ji-ki,i,  ^ji+ki+i,i)  ■  (24) 

The  support  extension  is  the  interior  of  the  set  formed  by  the  union  of  the  supports 
of  all  B-spline  basis  functions  whose  support  intersects  Q.  Note  that  each  element 
belongs  to  the  support  extension  of  at  most  n;=i^...^rf(2A;;  -|-  1)  elements. 


4.3  Piecewise  Smooth  Functions,  Geometrical  Mappings,  and 
Physical  Mesh  Entities 

On  the  parametric  mesh  we  define  the  space  of  piecewise  smooth  functions  with 
interelement  regularity  given  by  the  vectors  Oi, . . . ,  as 

=  (Mk) . 
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iOO 


(25) 


Precisely,  a  function  in  is  a  function  whose  restriction  to  an  element  Q  G  M.h 

admits  a  extension  in  the  closure  of  that  element  and  which  has  contin¬ 
uous  derivatives  with  respect  to  the  Ith  coordinate  along  the  internal  mesh  faces 
{{xi,...,Xd)  :  xi  =  7^  1]  for  alH«  =  2, . . . ,  m;  -  1  and 

ji'  =  1, . . .  —  1.  Note  immediately  that  any  function  lying  in  the  B-spline  space 

also  lies  in 

Unless  specihed  otherwise,  we  assume  throughout  the  rest  of  the  paper  that  the 
physical  domain  hi  C  can  be  exactly  parametrized  by  a  geometrical  mapping 
F  ;  n  — )■  n  belonging  to  o,^)  with  piecewise  smooth  inverse.  We  further 

assume  that  the  physical  domain  hi  is  simply  connected  with  connected  boundary 
dVt  and  the  geometrical  mapping  is  independent  of  the  mesh  family  index  h.  A 
geometrical  mapping  meeting  our  criteria  could  be  defined  utilizing  B-splines  or  Non- 
Uniform  Rational  B-Splines  (NURBS)  on  the  coarsest  mesh  Aiho-  For  examples  of 
such  mappings,  see  Chapter  2  of  [13] .  NURBS  mappings  are  especially  useful  as  they 
can  represent  many  geometries  of  scientihc  and  engineering  interest  and  are  the  main 
tools  employed  in  Computer  Aided  Design  (CAD)  software.  Later  in  this  paper,  we 
will  utilize  a  polar  mapping  to  solve  a  flow  problem  on  a  cylindrical  geometry.  The 
geometrical  mapping  F  naturally  induces  a  mesh 

lCh  =  {K-K=¥{Q),QeMh}  (26) 

on  the  physical  domain  D.  We  dehne  for  each  element  K  &  JCh  Si  physical  mesh  size 

=  ||-DF||2,oo(Q)hQ  (27) 

where  Q  is  the  pre-image  of  K,  and  we  also  dehne  the  support  extension  K  =  F((5). 
We  dehne  for  a  given  mesh  the  global  mesh  size 

h  =  max  {hx,  K  G  ICh}  ■ 

Note  that  as  the  parametric  mesh  family  {■Mh}h<ho  locally  quasi-uniform  and  the 
geometrical  mapping  F  is  independent  of  the  mesh  family  index  h,  the  physical  mesh 
family  is  also  locally  quasi-uniform.  We  refer  to  the  physical  domain  D  and 

its  pre-image  D  interchangeably  as  the  patch.  It  should  be  noted  that,  in  general,  the 
domain  D  cannot  be  represented  using  just  a  single  patch.  Instead,  multiple  patches 
must  be  employed.  We  will  discuss  further  the  multi-patch  setting  in  Section  7. 

We  dehne  on  the  parametric  mesh  a  set  of  mesh  faces  J^h  =  {F’}  where  F  is  a  face 
of  one  or  more  elements  in  M.h-  We  dehne  the  physical  set  of  mesh  faces  as 

:Fh  =  {F=Y{F)-.Feh] 

and  we  dehne  the  boundary  mesh  to  be 

Vh  =  {F  eFh-.F  dd^l}. 
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By  construction, 


dQ  =  Uperh-^- 

Note  that  for  each  face  F  ETh  there  is  a  unique  K  E  Kh  such  that  F  is  a  “face”  of 
K  (in  the  sense  that  F  is  the  image  of  a  face  of  Q,  the  pre-image  oi  K).  We  hence 
dehne  for  such  a  face  the  mesh  size 


hp  —  hp. 

One  may  also  dehne  hp  to  be  the  wall-normal  mesh-size  as  is  done  in  [3].  Such  a 
dehnition  is  more  appropriate  when  stretched  meshes  are  utilized  in  the  presence  of 
layers. 

Throughout  the  paper,  we  will  utilize  the  terminology  “a  constant  independent 
of  h” .  When  we  employ  such  terminology,  we  simply  indicate  that  the  constant 
will  not  depend  on  the  given  mesh  and,  in  particular,  its  size.  The  constant  may, 
however,  depend  on  the  domain,  the  shape  regularity  of  the  parametric  mesh  family, 
the  polynomial  degree  and  smoothness  of  the  employed  B-spline  spaces,  and  global, 
mesh-invariant  measures  of  the  parametric  mapping. 


5  Discretization  of  Velocity  and  Pressure  Fields 

In  this  section,  we  dehne  the  B-spline  spaces  which  we  will  utilize  to  discretize  the 
velocity  and  pressure  helds  appearing  in  the  steady  Navier-Stokes  problem.  These 
spaces  are  motivated  by  the  recent  theory  of  isogeometric  discrete  differential  forms 
[9,  10]  and  may  be  interpreted  as  smooth  generalizations  of  Raviart-Thomas  elements 
[32].  We  hrst  dehne  our  discrete  velocity  and  pressure  spaces  on  the  parametric 
domain  =  (0, 1)'^  and  then  dehne  discrete  spaces  on  the  physical  domain  using 
divergence-  and  integral-preserving  transformations.  We  hnish  this  section  with  a 
presentation  of  local  approximation  estimates  and  trace  inequalities  for  our  discrete 
velocity  and  pressure  spaces.  For  a  more  in-depth  discussion  of  the  discrete  velocity 
and  pressure  spaces  used  in  this  paper,  see  Section  5  of  [19]. 


5.1  Discrete  Spaces  on  the  Parametric  Domain 

Using  the  notation  of  the  previous  section  and  assuming  that 

a  :=  min{|Q:j|  :  /  =  1, . . . ,  d}  >  1, 

we  dehne  the  following  two  spaces: 

2, 
3, 


V„:= 


Ct  1 ,  Ct  2  —  1 , «  3  —  1 


*^Cti  ,Ct2  — 1 


Qfcl  — 1,/C2 
^  — 1,0:2 


X 


qki  —  l,k2,k3  —  l 
*^Oi  — 1,02  >0:3  — 1 


X 


qki  —  l^k2  —  liks 
*^01—1,02  —  1,03 


if  d  = 
if  d  = 
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if  d  =  2, 
if  d  =  3. 


qki  —  l,k2  —  l 
— 1,0:2  — 1 

qki  —  l^k2  —  l,k^  —  l 
*^01—1,02  —  1,03  —  1 


The  space  Vh  comprises  our  set  of  discrete  velocity  fields  while  Qh  comprises  our 
set  of  discrete  pressure  helds.  Note  that  as  a  >  1,  our  discrete  velocity  helds  are 
H ^-conforming.  If  we  allow  a  =  0,  our  spaces  collapse  to  standard  Raviart-Thomas 
mixed  finite  elements  [32],  In  order  to  deal  with  no-penetration  boundary  conditions, 
we  make  use  of  the  following  constrained  discrete  spaces: 


Vo,h  :=  <1  v/i  G  ;  v/j  ■  n  =  0  on 


Qo,h  ■=  Qh  -  /  g/idx  =  o 


in 


Above,  n  denotes  the  outward-facing  normal  to  dQ.  As  specihed  in  the  introduction, 
we  choose  to  enforce  no-slip  boundary  conditions  weakly  using  Nitsche’s  method  [30]. 
Due  to  the  special  relationship  given  by  (18),  the  spaces  Vo,/i  and  Qo,h  along  with  the 
parametric  divergence  operator  form  the  bounded  discrete  cochain  complex 


Vo,/.  Q, 


0,/i 


where  div  is  the  divergence  operator  on  the  unit  cube  D.  In  fact,  we  have  a  much 
stronger  result  due  to  the  results  of  [9]. 

Proposition  5.1.  There  exist  stable  projection  operators  If^  :  Ho(div;r2)  — )■  Vq,/. 
and  Ilg  ;  — )■  Qq,/.  such  that  the  following  diagram  commutes: 


Ho  (div;  Q) 


div 


V, 


div 


0,h 


n't 


.  0 


(28) 


0,h- 


Furthermore,  there  exists  a  positive  constant  Cu  independent  of  h  such  that 

l|np,yilHi(fi)  -  Vv  e  H„(div;S2)  n  H"(S)). 


(29) 


5.2  Discrete  Spaces  on  the  Physical  Domain 

To  dehne  our  discrete  velocity  and  pressure  spaces  on  the  physical  domain,  we  intro¬ 
duce  the  following  pullback  operators: 


t„(v)  :=  det(DF)(DF)“^(voF),  vGHo(div;D)  (30) 

ip{q)  :=  det  (DF)  (g  o  F) ,  g  G  Lo(^)  (31) 
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where  DF  is  the  Jacobian  matrix  of  the  parametric  mapping  F.  The  push-forward 
given  by  (30),  popularly  known  as  the  Piola  transform,  has  two  important  properties: 
(i)  it  preserves  the  nullity  of  normal  components,  (ii)  it  maps  divergences  to  diver¬ 
gences.  On  the  other  hand,  the  push-forward  given  by  (31)  has  the  property  that  it 
preserves  the  nullity  of  the  integral  operator.  Due  to  these  properties,  we  have  the 
following  commuting  diagram: 


Ho(div;  Q) 


Ho(div;  D) 


div 


(32) 


div 


This  motivates  the  use  of  the  following  discrete  velocity  and  pressure  spaces  in  the 
physical  domain: 

Vo,h  :=  |v  G  Ho(div;  D)  :  t„(v)  G  Vo,/i|  , 


Qo,h  ■—  ^ 

Furthermore,  we  dehne  the  projectors  :  Ho(div;r2)  — )■  Vo,h  and  — )■ 

Qo,/i  via  the  compositions 


Employing  the  preceding  results  and  terminology  as  well  as  the  smoothness  properties 
of  the  parametric  mapping  F ,  we  arrive  at  the  following  proposition. 

Proposition  5.2.  The  following  diagram  commutes: 


Ho(div;  fl) 


div 


TTO 


Qh 


(33) 


V, 


div 


0,/i 


^  Qo,h- 


Furthermore,  there  exists  a  positive  constant  Cu  independent  of  h  such  that 

l|nJ.,v||Hi,n)  <  C^lIvllHiin,,  Vv  e  H„(div;fi)  n  H'(fi). 


(34) 


We  immediately  have  an  inf-sup  condition  for  our  discrete  velocity /pressure  pair. 

Proposition  5.3.  There  exists  a  positive  constant  (3  independent  of  h  such  that  the 
following  holds:  for  every  qh  G  Qo,/i;  there  exists  a  ^  '^o,h  such  that: 

divvfe  =  qh  (35) 


12 


and 


(36) 


Hence, 

inf 

<?hSQo.(i 

Qh¥=0 

Proof.  See  the  proof  of  Proposition  5.3  in  [19].  □ 

We  also  have  the  following  result. 

Proposition  5.4.  If^h  G  Vo,/i  satisfies 

{dw\h,qh)LHn)  =  0,  G  Qo,/i,  (38) 

then  divv/j  =  0. 

Proof.  The  proof  holds  trivially  as  div  maps  Vo,/i  onto  Qo,/i-  n 

Hence,  by  choosing  Vo,/i  and  Qo,/i  as  discrete  velocity  and  pressure  spaces,  we 
arrive  at  a  discretization  that  automatically  returns  velocity  fields  that  are  pointwise 
divergence- free. 

5.3  Approximation  Results  and  Trace  Inequalities 

Let  us  dehne 

k'  =  min  \ki  —  1\ .  (39) 

1=1,.. .,d 

Note  that  the  discrete  velocity  and  pressure  spaces  Vo,/i  and  Qo,/i  consist  of  mapped 
piecewise  polynomials  which  are  complete  up  to  degree  k'.  Hence,  k'  may  be  thought 
of  as  the  polynomial  degree  of  our  discretization  technique.  The  following  result 
details  the  local  approximation  properties  of  our  discrete  spaces.  Its  proof  may  be 
found  in  [9]. 

Proposition  5.5.  Let  K  ^  JCh  and  K  denote  the  support  extension  of  K.  For 
d<j<s<k'  +  l,  we  have 

VveH‘(A')nH„(div;f!)  (40) 

\'}-K,.l\H.{K)<Ch‘ph\\„.iK),  Vg  e  44>(/0  n  L§(f!)  (41) 

where  C  denotes  a  positive  constant,  possibly  different  at  each  appearance,  independent 
ofh. 


(divvft,g/,)i2(n)  o 
sup  ^ ^ ^ ^ - >  fi. 

v/iSVo./i 
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Hence,  our  discrete  spaces  deliver  optimal  rates  of  convergence  from  an  approxi¬ 
mation  point  of  view.  We  will  also  need  the  following  trace  estimate  in  our  ensuing 
mathematical  analysis.  Its  proof  can  be  found  in  [16]. 

Proposition  5.6.  Let  K  E  ICh  and  Q  =  Then  we  have 

II  1 1  V/i  1 1  jjl(K),  Vv/j  G  Vo,h  (42) 

where  C trace  denotes  a  positive  constant  independent  of  h. 

In  [18],  it  was  shown  that  Proposition  5.6  holds  for  B-splines  and  parametric  hnite 
elements  with  Ctrace  ~  (^0^-  However,  our  numerical  experience  has  indicated  that 
a  corresponding  global  trace  inequality  holds  with  Ctrace  ~  k'  if  B-splines  of  maximal 
continuity  are  utilized.  This  allows  us  to  select  a  smaller  penalty  parameter  when 
employing  Nitsche’s  method.  As  we  will  see  in  the  next  section,  our  convergence 
estimates  scale  inversely  with  the  square  root  of  Nitsche’s  penalty  parameter.  Hence, 
we  want  to  select  Nitsche’s  penalty  parameter  as  small  as  possible. 


6  The  Discretized  Problem 

In  this  section,  we  approximate  the  homogeneous  steady  Navier-Stokes  problem  using 
the  discrete  velocity  and  pressure  spaces  introduced  in  the  previous  section.  We  prove 
continuity,  stability,  and  a  priori  estimates  for  our  discretization  scheme  in  the  single 
patch  setting  under  a  smallness  condition,  and  we  explicitly  track  all  of  our  estimates’ 
dependencies  on  the  viscosity  and  Nitsche’s  penalty  parameter. 


6.1  Variational  Formulation 


We  begin  this  section  by  presenting  a  discrete  variational  formulation  for  the  steady 
Navier-Stokes  problem.  Since  members  of  Vo,/i  do  not  satisfy  homogeneous  tangential 
Dirichlet  boundary  conditions,  we  resort  to  Nitsche’s  method  to  weakly  enforce  no-slip 
boundary  conditions.  Dehning  the  bilinear  form 


/Cft(w,  v)  =  A;(w,  v) 


2Z/  ((VW)n)  ■  w  -|-  ((V®w)  n) 


a 


pen 


hi 


w  ■  V 


ds 


(43) 

where  Open  >  1  is  a  chosen  positive  penalty  constant,  our  discrete  formulation  is 
written  as  follows. 


(G)<^ 


Find  Uh  G  Vo,/i  and  ph  G  Qo,h  such  that 

kh{uh,  v/i)  -t-  c{uh,  Uh,  Vft)  -  b{ph,  v/j)  -h  b{qh,  u^) 
for  all  Vfe  G  Vo,/i,  qh  e  Qo,/i- 


(fW/i)L2(o)  (44) 


We  have  the  following  lemma  detailing  the  consistency  of  our  numerical  method. 
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Lemma  6.1.  Suppose  that  (u,p)  is  a  solution  of  {W)  satisfying  the  regularity  con¬ 
dition  u  G  for  some  e  >  0.  Then: 

/c/i(u,  v/j)  +  c(u,  u;  Vft)  -  h{p,  \h)  +  b{qh,  u)  =  (f,  (45) 

for  all  e  Vo,h  and  qh  G  Qo.ft- 


Proof.  We  trivially  have 

b{qh,  u)  =  0,  Wqh  G  Qqa- 

Now  let  Vft  G  Vo,h-  By  the  Sobolev  trace  theorem,  the  assumption  u  G 
guarantees  that  (V^u)  n  is  well-dehned  along  dQ  and  (V®u)  n  G  {L‘^{dQ)Y.  Hence, 
the  quantity  kh{u,  \h)  is  well-dehned.  Utilizing  integration  by  parts  and  the  fact  that 
u  satishes  homogeneous  Dirichlet  boundary  conditions  and  v/^  satishes  homogeneous 
normal  Dirichlet  boundary  conditions,  we  have 


kh{u,  \h)  +  c(u,  u;  \h)  -  b{p,  v/^) 


(V  ■  (u  O  u) 


V  ■  (2z/V®u)  -h  gradp)  ■  v/^ 


This  completes  the  proof  of  the  lemma. 


□ 


We  have  the  following  relationship  between  the  exact  solution  and  a  numerical 
solution. 

Corollary  6.1.  Let  {uh,Ph)  denote  a  solution  of  {G),  and  let  (u,p)  denote  a  solution 
of  iW)  satisfying  the  regularity  condition  u  G  for  some  e  >  0.  Then: 

kh{u  -  u/i,  v/j)  -h  c(u,  u,  \h)  -  c(uh,  u/i,  \h) 

-b{p-ph,Vh)  +  b{qh,n-Uh)  =  0  (46) 


for  all  Vh  G  Vo,h  and  qh  G  Qo.ft- 

Finally,  by  Proposition  5.4,  we  have  the  following  lemma. 

Lemma  6.2.  Let  {uh,Ph)  denote  a  solution  of  (G).  Then: 

divuft  =  0  (47) 

Weak  imposition  of  no-slip  boundary  conditions  allows  our  methodology  to  default 
automatically  to  a  compatible  discretization  of  Euler  how  in  the  setting  of  vanishing 
viscosity.  Moreover,  for  large  Reynolds  number  hows,  there  is  a  sharp  boundary  layer 
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in  the  vicinity  of  walls.  Utilzing  Nitsche’s  method  allows  us  to  account  for  these  layers 
in  a  stable  and  consistent  manner  without  having  to  directly  resolve  them  [2,  3,  4],  In 
fact,  Nitsche’s  method  can  be  interpreted  as  a  variationally  consistent  wall  model.  To 
better  see  this  interpretation,  let  us  formally  rewrite  our  discrete  variational  equations 
as 


/T;V^Vftdx-^  /  Q-\hds  +  c{uh,Uh-,\h) -b{Ph,^h) +  b{qh,Uh)  = 

(48) 


where  T  is  a  symmetric  tensor  satisfying 


T  :  Wdx  = 


I  2z/V®Uft  :  Wdx  —  ^  j  2z/u/j  ■  (Wn)  ds 
Jn  pgP  J F 


=  /  2uuh  ■  divWdx  (in  the  sense  of  distributions)  (49) 


in 


for  symmetric  tensors  W  with  well-dehned  normal  trace  and  Q  a  vector  satisfying 


=  2z/  (V^u/i)  n  - 


''pen 

IF 


(50) 


Above,  T  is  a  weakly  dehned  viscous  stress  tensor  and  Q  is  the  resultant  viscous 
boundary  traction  vector.  In  the  event  that  the  no-slip  boundary  condition  is  met 
exactly,  we  recover  T  =  2z/V^U/i  and  Q  =  2z/  (V^u/j)  n.  Otherwise,  the  dehnitions 
of  T  and  Q  are  changed  accordingly.  The  tangential  component  of  Q  given  by  (50), 
denoted  Qtang,  and  calculated  as 


Oiona  =  Q  -  (Q  ■  n)  n, 


(51) 


is  the  effective  wall  shear  stress  vector. 

As  the  discrete  velocity  held  satishes  the  no-penetration  boundary  condition  strongly, 
the  vector  Q  is  equal  to  the  discrete  shear  stress  2z/  (V^U/^)  n  plus  an  additional  wall 
shear  stress  term  Q"*"  in  the  direction  tangent  to  the  wall.  Specihcally,  we  have 


where 


=  —u 


*2 


u*^  = 


''pen 


h 


Uh 

1  II 

(52) 

|Uh|| 

l|Uh|| 

(53) 

For  under-resolved  how  simulations,  the  magnitude  of  (V^U/^)  n  in  the  direction  tan¬ 
gent  to  the  wall  is  relatively  small  and,  as  such,  the  tangential  component  of  Q  is 
dominated  by  Q"*".  In  this  sense,  becomes  a  model  for  the  wall  shear  stress.  As 
the  mesh  is  rehned  and  the  how  is  resolved,  Q"*"  — )■  0.  The  above  interpretation  allows 
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US  to  design  physically  motivated  penalty  values  for  Nitsche’s  penalty  parameter.  No¬ 
tably,  u*  may  be  interpreted  as  the  friction  velocity.  By  specifying  the  value  of  u* 
using  Spalding’s  law  of  the  wall  [36],  we  recover  a  standard  wall  model  for  under¬ 
resolved  flow  simulations.  For  more  on  this  approach,  see  Section  3  of  [3].  We  recall, 
however,  that  the  numerically  inspired  (53)  produced  results  of  the  same  quality  as 
the  u*  given  by  Spalding’s  physically  inspired  law  of  the  wall. 


Remark  6.1.  If  we  wish  to  impose  non-homogeneous  tangential  Diriehlet  (e.g.,  pre¬ 
sented  slip)  boundary  conditions,  we  must  add  the  following  expression  to  the  right 
hand  side  of  our  discrete  formulation: 


/v(vh)  = 


■ 


2iy  (  -  ((VW/,)  n)  ■  ubc  +  ds  (54) 


where  vlbc  is  a  prescribed  vector  function  defined  on  dQ.  If  we  also  wish  to  impose 
non-homogeneous  normal  Diriehlet  (e.g.,  prescribed  penetration)  boundary  conditions, 
we  must  impose  these  strongly  and  add  the  following  expression  to  the  left  hand  side 
of  our  discrete  formulation: 


Cf/vy(u/i,  V/j)  — 


E 

FGFh 


{ubc  ■  n) ,  U/J  ■  \hds 


(55) 


and  the  following  expression  to  the  right  hand  side  of  our  discrete  formulation: 


fuwi'^h)  —  —  /  (usc  ■  n)_  Ubc  '  ^hds 

Fer^ 


(56) 


where 


and 


(use  ■  n) ,  = 


Ubc  ■  n  if  u^c  ■  n  >  0 

0  otherwise 


{ubc  ■  n)_  = 

These  additional  terms  correspond  to  upwinding. 


Ubc  ■  n  if  ubc  ■  n  <  0 

0  otherwise. 


6.2  Well-Posedness  for  Small  Data 

We  now  prove  that  our  discrete  formulation  is  well-posed  under  a  smallness  condition. 
Onr  method  of  proof  mimics  that  of  the  continnous  problem  (see  Theorem  10.1.1  of 
[22]).  To  begin,  let  us  dehne  the  following  mesh-dependent  norm: 

l|v|lh  :=  |v|h1(o)  +  hpW  (VW)n||2^2(^)),  +  ^  (57) 

Fer,  Fer, 


17 


We  will  also  denote 


Vo,/i  :=  {v/i  G  Vo, ft  :  divvft  =  0}  =  {vft  G  Vo, ft  ;  &(gft,  Vft)  =  0,  Vgft  G  Qo,ft}  • 


To  proceed,  we  need  to  call  upon  stability  and  continuity  results  that  were  proven  in 
[19]  for  divergence-conforming  B-spline  discretizations  of  the  Darcy-Stokes-Brinkman 
equations.  These  results  hinge  upon  two  assumptions  regarding  the  size  of  Open-  First, 
in  light  of  Proposition  5.6,  we  assume  that 


^pen  Poin^ Korn 


V/,1 


ViF  G  /Cft,  Vft  G  V 


o,ft 


(58) 


where  Cpoin  is  the  Poincare  constant  appearing  in  (10)  and  Ckotu  is  the  positive 
constant  associated  with  the  following  Korn’s  inequality  [6]: 


|w|Ll(0)  <  CKorn  (ll  V"w ||  II W ||  ,  Vw  G  H^(fl). 

Second,  we  assume  that 

C'pen>4ho|9fi|-'/(‘'-l)  (59) 

where  ho  is  the  mesh  size  of  the  coarsest  mesh  /Co  and  |(9r2|  denotes  the  length  of  dQ  for 
d  =  2  and  the  area  of  dfl  for  d  =  3.  This  second  assumption  is  necessary  as  rotation 
modes  carry  zero  energy.  Hence,  weak  boundary  conditions  are  needed  to  control  these 
modes  in  rotationally  symmetric  (or  near  rotationally  symmetric)  conhgurations.  As 
such  conhgurations  are  of  signihcant  engineering  interest,  we  believe  that  any  analysis 
results  should  cover  these  situations.  Note  that  a  constant  Open  satisfying  the  above 
assumption  need  not  depend  on  h  or  u.  Rather,  it  only  needs  to  depend  on  the 
size  of  the  domain,  the  polynomial  degree  and  smoothness  of  the  discretization,  the 
parametric  shape  regularity,  and  global,  mesh-invariant  measures  of  the  parametric 
mapping. 


Corollary  6.2. 

b{p,  v)  < 


4(Wft,  Wft)  > 


Assume  (58)  and  (59)  are  satisfied.  Then  we  have 


2z/Cco„t||w||ft||v||ft, 

lblU2(o)||v||ft 

2z/Ccoerc||Wft||^, 


Vw,  V  G  Vo,ft  ©  (Hj(fl)  n  HV2+^(fl))  (60) 

Vp  G  T^(H),v  G  Vo,ft  ©  (Hj(fl)  nHV2+^(fl)) 

(61) 

Vw/i  G  Vo, ft  (62) 


where  e  >  0  is  an  arbitrary  positive  number  and  Ccont  and  Ccoerc  are  positive  constants 
which  are  independent  of  h,  o,  Open,  and  e.  Furthermore,  we  have 


inf  sup 

,h  V/iEVq  h 


(divvfe,gfe) 

||Vft||v(ft)|kft||Q 


>/3. 


(63) 


where  fi  is  a  positive  constant  independent  of  h  and  v  which  asymptotically  scales 
inversely  with  the  square  root  of  Open- 
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We  will  also  need  the  following  two  lemmas.  The  first  lemma  gives  a  Lipschitz 
continuity  result  for  the  trilinear  form  c(-,  •;  ■).  This  result  hinges  upon  Sobolev  em¬ 
beddings  which  exist  because  our  domain  is  Lipschitz.  The  second  lemma  gives  a 
semi-coercivity  result  for  the  trilinear  form  c(-,  •;  •). 

Lemma  6.3.  There  exists  a  constant  Cq  only  dependent  on  hi  such  that 

c(wi,x;v)  -c(w2,x;v)  <  (TqIwi  -  W2|h1(o)  |x|h1(o)  IvIhTJ^)  (64) 

for  all  wi,  W2,  X,  V  G  H^(r2)  fl  Ho(div;  hi) . 


Proof.  Let  wi,W2,x, v  G  H^(r2)  fl  Ho(div;r2)  be  arbitrary.  Note  that  since  dil  is 
Lipschitz,  we  have  the  continuous  embedding 

H\n)  -G  L\n). 

By  linearity  and  the  Cauchy- Schwarz  inequality,  we  can  then  write 
c(wi,  x;  v)  -  c(w2,  x;  v)  =  -  ((wi  -  W2)  ®  x,  Vv)(^2(n))dxd 

<  ||Wi  —  W2||L4(f^)||x||L4(Q)||  Vv||(i2(f2))dxd. 

Let  Cembed  deuote  the  positive  embedding  constant  dependent  only  on  the  domain  hi 
such  that 

l|y||L4(0)  <  Cemhed\\y\\'H^(n)-i  Vy  G  H  (hi). 

Then  we  have 

c(wi,x;v)  -c(w2,x;v)  <  -  W2||h1(0)  ||x||h1(0)  ||  Vv||(i2(f2))dxd. 

The  lemma  follows  with  Cq  =  where  Cpoin  is  the  Poincare  constant  ap¬ 
pearing  in  (10).  □ 


Lemma  6.4.  Suppose  w,v  G  H^(r2)  fl  Ho(div;r2)  such  that  divw  =  0.  Then 

c(w,  v;  v)  =  0.  (65) 


Proof.  We  write 

Since  divw  =  0,  we  have 


(w  ®  v)  :  Vvdx. 


div  (w|vp)  dx. 


The  lemma  is  then  simply  a  result  of  the  divergence  theorem. 


□ 
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With  Corollary  6.2  and  Lemmata  6.3  and  6.4,  we  can  prove  the  following  propo¬ 
sition  which  gives  the  well-posedness  of  the  linearized  Oseen  problem. 

Proposition  6.1.  Let  w  G  H^(r2)  fl  Ho(div;r2)  such  that  divw  =  0,  and  assume 
(58)  and  (59)  are  satisfied.  Then  the  following  problem  has  a  unique  solution:  find 
(xft,  rh)  e  Vo,h  X  Qo,h  such  that 


khi^h,  Vft)  c(w,  Xft;  v/i)  -  b{rh,  v^)  -h  b{qh,  x/^)  =  (f,  Vfe)L2(f2)  (66) 

for  all  \h  £  Vo,/i,  qh  G  Qo,h-  Furthermore,  the  unique  solution  of  the  linearized  Oseen 
problem  satisfies  divx/^  =  0  and  the  bound 


L2(Q) 


(67) 


where  Cpoin  is  the  Poincare  constant  appearing  in  (10)  and  Ccoerc  is  the  coercivity 
constant  appearing  in  Corollary  6.2. 


Proof.  Existence  and  uniqueness  are  a  direct  result  of  Brezzi’s  theorem,  Corollary  6.2, 
and  Lemmata  6.3  and  6.4.  Since  divVo,/i  =  Qo,h,  we  automatically  have  divx^  =  0. 
To  prove  the  a  priori  bound,  we  write  using  (66)  and  the  coercivity  of  kh{-,  •) 


2iy  Ccoerc 
1 


2izCc 


kh(Kh,Xh) 

((f,X/.)L  2(0)  -  c(w,  X/,;  x/,)  -  b{rh,  x^)) 


Using  Lemma  6.4  to  set  c(w,x;i;x/i)  =  0  and  divx/^  =  0  to  set  h{r-h,^i() 
complete  the  proof: 


xkWI  < 
< 

< 


2vC coerc 
1 

2vC coerc 

Cpoin 
2vC cnerc 


(f)  X/i)l2(q) 

l|f||L2(0)  I|x/i||l2(0) 
l|f||L2(0)  |Xh|Hi(!^)- 


0,  we  can 


□ 


We  are  now  ready  to  establish  well-posedness  results  for  the  full  Navier-Stokes 
problem.  To  do  so,  we  will  attempt  to  obtain  the  solution  to  the  discrete  Navier-Stokes 
problem  through  the  iterative  solution  of  a  sequence  of  Oseen  problems.  Notably, 
given  uq  G  Vq  ^,  we  seek  the  limit  of  the  following  iterative  solution  process:  for 
i  =  1,  2,  3, . . .  hnd  (uj,pj)  G  Vo,/i  x  Qo,/i  such  that 

kh{vLi,  Vfe)  4-  c(uj_i,  Ui,  \h)  -  b{pi,  Vft)  -4  b{qh,  u*)  =  (f,  v;i)L2(f2)  (68) 
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for  all  Vft  G  Vo,/i,  Qh  G  Qo,h-  Well-posedness  hinges  upon  the  convergence  of  this 
iterative  solution  process,  and  convergence  of  this  process  is  further  contingent  upon 
a  small  data  constraint. 


Theorem  6.1.  Assume  (58)  and  (59)  are  satisfied,  and  further  assume  that 


CaC 


O^Poin 


|f|lL2(0)  <  1 


(69) 


where  Cq  is  the  eontinuity  constant  appearing  in  Lemma  6.3,  Cpoin  is  the  Poincare 
constant  appearing  in  (10),  and  Ccoerc  is  the  coercivity  constant  appearing  in  Corollary 
6.2.  Then  Problem  (G)  has  a  unigue  solution  {uh,Ph)  G  Vq  ,h  X  Qo,h-  Furthermore, 
the  unigue  solution  satisfies 


IlUftlU  < 

lb/i||L2(0) 


c 


Poin 


2vC. 

<r' 


coerc 

c^c 


Il2(o) 


0^  Poin 


+ 


a 


cont 


(2z/)2CL.e  a 


+  l]c_ 


-'Poin 


lL2(f^) 


(70) 

(71) 


where  fi  and  Ccont  are  the  inf-sup  and  continuity  constants  appearing  in  Corollary 

6.2. 


Proof.  Begin  by  defining  S  :  Vo,/i  — )■  Vo,h  to  be  the  nonlinear  operator  which  returns 
the  divergence-free  velocity  solution  of  (66)  given  a  divergence-free  velocity  held  G 
Vq  /j.  Note  that  by  Proposition  6.1 


l|S(w.)lk<;5%^^l|f| 


Lho)- 


Therefore,  the  nonlinear  operator  S  maps  Vg  ^  into 


Bh  :=  <  w/i  G  Voh  ■  llwftlU  < 


C 


Poin 


2vCr. 


Il2(o) 


Now,  let  wi,W2  G  and  Wi  =  S'(wi),  W2  =  S'(w2).  By  using  the  coercivity  of 
kh{-,-)  given  by  Corollary  6.2  we  have 

2z/Ccoerc||wi  -  W2||^  <  hifiVi  -  W2,  Wi  -  W2).  (72) 


By  (66),  we  have 


kh(wi,  Wi  -  W2)  =  -C(wi,  Wi,  Wi  -  W2)  -h  (f,  Wi  -  W2)l2(0) 


and 


kh(w2,  Wi  -  W2)  =  -c(w2,  W2,  Wi  -  W2)  -h  (f,  Wi  -  W2)l2(0) 
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since  div  (wi  —  W2)  =  0.  Hence, 

2z/Ccoerc||wi  -  W2\\l  <  c(w2,  W2,  Wi  -  W2)  -  c(wi,  Wi,  Wi  -  W2) 

=  -C(W2,  Wi  -  W2,  Wi  -  W2) 

+  C(W2,  Wi,  Wi  -  W2)  -  C(wi,  Wi,  Wi  -  W2). 

By  the  continuity  and  coercivity  properties  given  in  Lemmata  6.3  and  6.4,  we  can 
write 


2z/Ccoerc||wi  -  W2||^  <  Co||w2  -  Wi  ||  ft  ||  Wi  || /^  ||  Wi  -  W2||ft. 

Since  wi  G  Hft,  we  thus  have 


Wi  -  W2||ft  <  fj,\\w2  -  Willft 


where  /x  is  precisely 


CqC Poin  IIP  II 

By  assumption  /i  <  1,  and  we  therefore  have  proved  S'  is  a  contractive  map.  By  the 
Banach  hxed  point  theorem,  the  nonlinear  problem 


Uft  =  S(uft) 


has  a  unique  solution  which  lies  in  Bh  and  is  precisely  the  discrete  velocity  solution 
(now  proven  unique)  of  Problem  {G).  Given  Uft,  uniqueness  of  the  discrete  pressure 
solution  ph  is  a  direct  result  of  the  inf-sup  condition. 

We  now  prove  the  stability  bounds.  The  bound  for  Uft  is  straight-forward  since 
Uft  G  Bh-  To  prove  the  bound  for  pft,  we  utilize  the  inf-sup  condition,  the  continuity 
of  kh{-,-),  the  continuity  of  c(-,  •;  •),  and  Poincare’s  inequality: 


P\\Ph\\L'^{n) 


<  sup 

VhSVo,;, 


=  sup 

vjiSVo,;! 


KPh,  Vfe) 
hhWh 

-(f,  Vft)L2(j2)  +  kh{uh,  Vft)  -h  c(uft,  Uft;  Vft) 

ilvhIU 


< 


l|f|lLhO)l|Vft||Lhn)  +2z/C'cont||Uft||fti|Vft||ft  +  C'o||Uft||2||vft||ft 

sup  - ^ ^ - 

VhSVo.h  l|V/i||/i 


< 


Cpom||f||L2(n)l|v/i|U  +  2z/Gcont||Uft||ft||Vft||ft  Gq || Uft || ^ || Vft || ft 
sup  - ^ ^ - 

VhSVo.h  W'^h]  h 


The  result  then  follows  by  invoking  the  bound  for  Uft. 


□ 
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To  the  best  of  our  knowledge,  the  above  theorem  is  the  first  of  its  kind  for 
quadrilateral-  and  hexahedral-based  spline  hnite  element  discretizations  written  in 
divergence-form.  Such  a  result  also  hypothetically  exists  for  Scott- Vogelius  discretiza¬ 
tions  on  tetrahedra  [41],  but  these  discretizations  are  limited  to  an  extremely  restric¬ 
tive  class  of  macro-element  meshes.  Well-posedness  gives  us  the  confidence  that,  at 
the  very  least,  our  formulation  gives  a  unique  solution  under  a  smallness  condition 
not  unlike  that  for  the  continuous  problem.  In  the  next  section,  we  will  show  that 
our  discrete  solution  converges  to  the  exact  solution  under  a  slightly  more  restrictive 
smallness  condition. 

It  is  interesting  to  note  that  the  proof  of  Theorem  6.1  guarantees  that  the  hxed 
point  iteration  given  by  (68)  converges  to  the  exact  solntion  for  any  initial  divergence- 
free  velocity  held.  Furthermore,  the  iterates  exhibit  the  linear  convergence  rate 

||u„+i  -  u^ll/i  <  /i”||ui  -  UoWh 


for 


= 


C'nC', 


0^  Poin 


|f||L2(0)  <  1- 


This  hxed-point  scheme  can  be  accelerated  using  a  Newton-Raphson  procedure. 


6.3  A  Priori  Error  Estimates  for  Small  Data 


We  are  now  ready  to  show  that  onr  discrete  solution  helds  converge  to  the  exact 
solution  helds  nnder  smallness  and  regnlarity  conditions.  Our  method  of  proof  largely 
mimics  that  of  Theorem  4.8  in  [11]  for  the  two-dimensional  Navier-Stokes  equations. 
However,  our  proof  is  more  straight-forward,  primarily  due  to  four  facts:  (1)  we 
employ  natively  divergence-free  discretizations,  (2)  we  employ  smooth  approximation 
spaces,  (3)  we  have  a  simpler  treatment  of  the  convection  operator,  and  (4)  we  do 
not  include  stress  as  an  auxiliary  variable. 

Our  hrst  error  estimate  reads  as  follows.  Note  that  it  is  explicit  in  the  mesh-size 
h,  the  dihnsivity  z/,  and  the  penalty  parameter  Open- 

Theorem  6.2.  Assume  (58)  and  (59)  are  satisfied,  and  further  assume  that 


max 


OnO, 


O^Poin 


OoO, 


^^Poin 


OnO 


’  Z/20e, 


O^Poin 


l|f|lL2(0)  <  1 


(73) 


where  Cq  is  the  eontinuity  eonstant  appearing  in  Lemma  6.3,  Cpoin  is  the  Poineare 
constant  appearing  in  (10),  Cq  is  the  domain- dependent  constant  appearing  in  The¬ 
orem  3.1,  and  Ccoerc  is  the  coercivity  constant  appearing  in  Corollary  6.2.  Let  (u,p) 
and  {uh,Ph)  denote  the  solutions  of  Problems  {W)  and  (C)  respectively.  Linder  the 
assumption  that  u  G  for  some  e  >  0,  we  have 


u-Uh\\h<{l  +  ‘2nj)  inf  ||u-Vfe||,j 

VhSVo,?! 


(74) 
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and 


lb-P/^||L2(^7)  <  (  1  + i  )  inf  \\p-qh\\LHn)  +  1\\u-Uh\\h  (75) 

\  jj  J  qh&Qo,h  jj 


where 


K  ~\~  1  “1“ 


1  \  CaC 


2Cr^ 


0^  Poin 
V 


7  = 


l|f||L2(Q)  <  2z/  {Ccont  +  2Ccoerc  +  C'coerc)  ,  (76) 
1 


2vCr. 


(77) 


and  fd  and  Ccont  (ire  the  inf-sup  and  continuity  constants  appearing  in  Corollary  6.2. 


Proof.  We  begin  by  proving  the  estimate  for  the  velocity  error.  Let  \h  ^  '^o,h  such 
that  divv/i  =  0.  By  the  coercivity  of  kh{-,  ■)  over  Vo,/n  we  have 

2z/Ccoerc||Vft  -  \lh\\l  <  kh{^h  “  U/j,  Vfe  -  Uh) . 

By  the  consistency  given  by  Corollary  6.1  and  the  divergence-free  condition  div(v/j  — 
u/j)  =  0,  we  have 

khiyh  -  u/i,  \h  -  Uh)  =  khiyh  -u,\h-  Uh) 

-  c(u,  u;  Vft  -  u,,)  c{uh,  u,,;  v,,  -  Uh). 

Then,  by  using  the  continuity  of  kh{-,-)j  we  can  write 

2uCcoerc\\^h  “  W||^  <  4(Vft  -  U,  Vft  -  U,,) 

-  c(u,  u;  \h  -  Uh)  +  c{uh,  Uh,  -  Uh) 

<  2uCcont\Wh  -  u|U||v,j  -  Uftll,,  -h  T  (78) 

where 


T  :=  -c(u,  u;  v/,  -  u,*)  c{uh,  u^;  \h  -  w). 
We  now  utilize  a  splitting.  Let  us  decompose 

T  =  Ti  +  Ts  +  Ts  +  T4 


where 


(79) 


Ti  =  -c(vfe,  u;  Vft  -  Uh)  +  c{uh,  u;  v/^  -  u^) 
T2  =  -c{uh,  v/i  -  u/,;  Vft  -  u,,) 

T3  =  c(yh,  u;  Vfe  -  Uh)  -  c(u,  u;  \h  -  u^) 

T4  =  c{uh,  ^h  -  u;  V/I  -  Uh). 
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By  Lemma  6.3,  we  can  write 


Ti  <  C'o|u|Hi(n)||vfe-u,,||^ 

T3  <  Co|u|Hi(n)||vfe  -  u/i||ft||v/i  -  u\\h 

T4  <  Co||u,,||,i||v/i  -  UfelUilVfe  -  u\\h 


and  by  Lemma  6.4,  we  have 


To  =  0. 


Theorems  3.1  and  6.1  immediately  give  the  bounds 


T3< 


V 

C^Cpoin  I 


lL2(0)l|V/i-Ufe|U||Vfe-U||;, 


V 

T4  <  ^°^^°'''||ffe||L2(0)||Vfe  -  Ufe||fe||Vfe  -  U||fe 
ZZ/Ocoer’c 

Substituting  (79),  (80),  (81),  (82),  and  (83)  into  (78),  we  obtain 

(1  -  a)  ||v/i  -  Uh\\h  <  «7l|vft  -  u\\h 


where 


C'nC' 


a  = 


O^Pom 


2Z/2C', 


coerc 

1 


Il2(o)) 


7  = 


2uCr. 


and 


K  2vC T  1  T 


2a, 


C^Cpoin  I 


Il2(o)- 


By  assumption,  a  <  1/2  and  hence 

||vfe  -  UhWh  <  2k7||v,i  -  u\\h. 


(80) 


(81) 

(82) 

(83) 


To  hnish  the  proof  for  the  velocity  error,  we  perform  a  sum  decomposition  as  follows: 
||u-u,i||/,<  inf  (||u  -  v/illfe  +  ||vft  -  Ufell/j) 

VhSVo.Ji 

<(1  +  2^7)  inf  ||u-v/i||,,. 


We  now  proceed  to  the  estimate  for  the  pressure  error.  Let  qh  G  Qo,/i-  By  the 
inf-sup  condition,  we  have 


hh-Vh\\L\n)<^  sup 

[j  VhGVo^h 


K^h  -Ph,^h) 

llvhik 


(84) 
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By  the  consistency  given  by  Corollary  6.1,  we  have,  for  any  v/^  G  Vo,/i, 

K^h  -  Phj  v/i)  =  b{qh  -  p,  Vfe)  +  kh{u  -  u/i,  \h) 

+  c(u,  u;  v/i)  -  c{uh,  Uh,  Vfe) 

=  b{qh  -  p,  Vfe)  +  kh{u  -  u/i,  \h) 

+  c(u,  u  -  u,,;  Vft)  +  c(u  -  U/I,  u^;  v/,). 

Our  continuity  estimates  from  Corollary  6.2  and  Lemma  6.3  then  give  the  bound 

K<lh  -  Ph,  V/i)  < 

{Uh  -  P||l2(0)  +  {‘2l^Ccont  +  Co  (|u|Hl(n)  +  ||Uh|U))  ||u  -  Uh\\h)  ||v/i|U. 

Theorems  3.1  and  6.1  give 

‘^^Ccont  +  Co  (|u|fjl(Q)  +  ||u/i||/i)  <  K, 

so  we  can  write 


K(!h  -  Ph,  Vfe)  <  {\\qh-  p||l2(o)  +  «||u  -  u/illft)  ilvftll,,. 

Substituting  the  above  expression  into  (84),  we  acquire  the  estimate 

1  Av 

\\qh-Ph\\L^n)  <  ^||gh-p||L2(n)  +  ^||u-u/,||,,. 

fj  fJ 

To  hnish  the  proof  for  the  pressure  error,  we  again  perform  a  sum  decomposition  as 
follows: 


||p-P/i||l2(o)  <  inf  {\\p- qh\\L^n)  +  \\qh-Ph\\L^n)) 

Qh^Qo,h 

+  \\p-qh\\L\n)+ '^\\u-Uh\ 

\  jj  J  qh&Qo,h  fj 


□ 

Our  next  error  estimate  gives  us  a  priori  convergence  estimates  that  are  optimal 
for  the  discrete  velocity  held  and  suboptimal,  by  one  order,  for  the  discrete  pressure 
held. 

Theorem  6.3.  Let  the  assumptions  of  Theorem  6.2  hold  true.  Furthermore,  let  (u,p) 
and  {uh,Ph)  denote  the  solutions  of  Problems  {W)  and  {C),  respectively.  Under  the 
assumption  that  u  G  and  p  G  for  some  j  >  1/2,  we  have 

||u  -  u/jII/j  <  Cu  (1  +  2^7)  h^||u||Ho+i(n)  (85) 
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and 


\\p-Ph\\LHn)  <  Cp 


K  , 


hmpWHHn)  +  ^||u  - 

P 


(86) 


for  s  =  min{A;',  j}  where  k  is  defined  by  (76),  7  is  defined  by  (77),  fi  is  the  diserete  inf- 
sup  constant,  Cu  is  a  positive  constant  independent  of  h  and  v  which  asymptotically 
scales  with  the  square  root  of  Open,  and  Cp  is  a  positive  constant  independent  of  h,  u, 


Proof.  We  first  prove  (85).  Recall  the  error  estimate  given  by  (74): 

||u  -  Ufe||/j  <  (1  +  2^7)  inf  ||u-v/,||/i. 

vjieVo./i 

Noting  divlly^u  =  divu  =  0,  we  can  choose  v/^  =  hly^u  in  the  above  expression 
to  obtain 

||u  —  u/ill/j  <  Cco \/Ti  T2  T3  (87) 

where  we  have  assigned  Cco  =  (1  +  2/1:7)  and 

Ti  =  |u- (88) 
T,  =  ^  AfII  (V*  (u  -  n».u))  n||y,y„.  (89) 

Ts  =  Cpenhp  i|u—  (90) 

F&Th 

To  handle  the  face  integral  in  (89),  we  recruit  the  multiplicative  trace  inequality 
for  fractional  Sobolev  spaces  [38]  and  Young’s  inequality  element-wise  to  obtain  the 
bound 

^  hpW  (V"  (u  -  n°^u))  n||2^2(^))d  < 

F&Vh 

iCtrcp)  ^  (\U  —  hly^uljj!^^^  -I-  h^|u  —  hly^  u| j 

K&Kh 

where  1/2  <  g  <  s  and  Ctrc,i  is  a  positive  constant  independent  of  h,  u,  and  Cpen-  To 
handle  the  face  integral  in  (90),  we  recruit  the  standard  continuous  trace  inequality 
element-wise  to  obtain  the  bound 


Cpenhp^\\u  ny^u||^^2('j7’))d  Y 


{Ctrc,2f  (Vl|u  -  n°^u||^2(^)  +  |u  -  UZ_  U 


2 

Vh  ‘^lHi(K) 


KeiCh 
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where  Ctrc,2  is  a  positive  constant  independent  of  h  and  u  which  varies  linearly  with 
the  sqnare  root  of  Open-  It  should  be  noted  the  two  constants  Ctrc,i  and  Ctrc,2  neces¬ 
sarily  depend  on  the  shape  regularity  of  the  mesh  family  parametric 

mapping  which  together  give  the  shape  regularity  of  the  mesh  family  See 

[18]  for  more  details.  Inserting  the  above  two  inequalities  into  (87)  and  then  applying 
Proposition  5.5,  we  immediately  acquire  the  bound 

||u  - 

for  Cu  a  positive  constant  independent  of  h  and  u  with  the  same  functional  depen¬ 
dency  on  the  penalty  parameter  as  Ctrc,2- 

The  proof  for  (86)  is  much  more  immediate.  Choosing  qh  =  hlg^p  in  the  error 
estimate  given  by  (75),  one  obtains 

lb  -  PhWm^)  <  (^1  +  lb  -  n^^p|b2(j^)  +  |lb  -  u,,|b  (9i) 

Inequality  (86)  follows  by  an  application  of  Proposition  5.5  to  bound  the  pressure 
interpolation  error.  □ 

Again,  to  the  best  of  our  knowledge.  Theorems  6.2  and  6.3  are  the  first  of  their 
kind  for  quadrilateral-  and  hexahedral-based  spline  hnite  element  discretizations 
written  in  divergence-form.  Note  that  we  have  obtained  error  estimates  which  are 
optimal  for  the  velocity  field  and  suboptimal,  by  one  order,  for  the  pressure  field 
under  a  smallness  condition  not  unlike  that  of  the  continuous  problem.  In  Section 
8,  we  will  employ  a  selection  of  problems  with  known  analytical  solutions  to  conhrm 
our  theoretical  convergence  rates.  Our  numerical  results  will  suggest  our  derived 
pressure  error  estimates  may  be  conservative.  The  analysis  presented  here  also  covers 
singular  solutions  typically  encountered  in  practice.  We  also  numerically  study  the 
effectiveness  of  our  method  for  a  singular  test  problem,  lid-driven  cavity  flow. 

7  Extension  to  Multi-Patch  Domains 

As  was  mentioned  previously,  most  geometries  of  scientihc  and  engineering  interest 
cannot  be  represented  by  a  single  patch.  Instead,  the  multi-patch  concept  must  be 
invoked.  We  assume  that  there  exist  Up  sufficiently  smooth  parametric  mappings 
Fj  ;  (0,  l)'^  — )■  such  that  the  subdomains 

1^2  F,:  ^  1 5  *  *  *  : 

n  = 
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are  non-overlapping  and 


Figure  1:  Example  multi-patch  construction  in 


We  refer  to  each  subdomain  (and  its  pre-image)  as  a  patch.  For  a  visual  depiction 
of  a  multi-patch  construction  in  see  Figure  1.  We  build  discrete  velocity  and 

pressure  spaces  over  each  patch  i  =  in  the  same  manner  as  in  the 

previous  sections  except  that  we  do  not  yet  enforce  boundary  conditions,  and  we 
denote  these  spaces  as  Vh(fii)  and  Qh{^i)- 

To  proceed  further,  we  must  make  some  assumptions.  First  of  all,  we  assume  that 
if  two  disjoint  patches  flj  and  Qj  have  the  property  that  dQi  fl  dQj  7^  0,  then  this 
intersection  consists  strictly  of  patch  faces,  edges,  and  corners.  More  succinctly,  two 
patches  cannot  intersect  along  an  isolated  portion  of  a  face  (or  edge)  interior.  Second, 
we  assume  that  the  mappings  are  compatible  in  the  following  sense:  if  two 

patches  flj  and  Qj  share  a  face,  then  Fj  and  Fj  parametrize  that  face  identically  up 
to  changes  in  orientation.  Third,  we  assume  that  if  two  patches  flj  and  Qj  share  a 
face,  the  B-spline  meshes  associated  with  the  patches  are  identical  along  that  face. 
This  guarantees  our  mesh  is  conforming.  Finally,  we  assume  for  simplicity  that  ki  = 
. . .  =  kd  =  k*  for  all  patches.  The  mixed  polynomial  degree  case  introduces  additional 
complications  that  are  beyond  the  scope  of  this  work.  We  would  like  to  note  that  all 
four  assumptions  hold  if  we  employ  a  conforming  NURBS  multi-patch  construction. 
See,  for  example.  Chapter  2  of  [13]. 
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We  define  our  global  discrete  velocity  and  pressure  spaces  as  follows: 

Vo,h  :=  {v;,  G  Ho  (div;  ff)  ;  G  Vi  =  1, . . . ,  np}  ,  (92) 

Qo,h  ■=  {(Ih  £  Lq  (ff)  :  ghlfii  G  Qh{fii),  Vi  =  1, . . .  ,np}  .  (93) 

The  space  Vo,/i  is  easily  constructed  due  to  our  preceding  four  assumptions  and  use 
of  open  knot  vectors.  Specifically,  we  set  to  zero  the  coefficient  of  any  basis  function 
whose  normal  component  is  nonzero  along  dQ,  and  along  shared  faces  between  patches 
we  (i)  equivalence  the  coefficients  of  any  basis  functions  whose  normal  components 
are  nonzero  and  equal  in  magnitude  and  direction  and  (ii)  set  the  coefficient  of  one  to 
minus  the  coefficient  of  another  for  any  basis  functions  whose  normal  components  are 
nonzero,  equal  in  magnitude,  and  opposite  in  direction.  We  note  that  this  is  precisely 
the  same  procedure  as  is  used  to  construct  Raviart-Thomas  spaces  on  conforming 
hnite  element  meshes.  We  simply  have  patches  instead  of  elements.  It  is  easily  shown 
that  the  spaces  Vo,/i  and  Qo,h,  along  with  the  divergence  operator,  form  the  bounded 
discrete  cochain  complex 

yo,h  - > 

However,  functions  in  Vo,/i  do  not  necessarily  he  in  H^(r2)  as  tangential  continuity 
is  not  enforced  across  patch  interfaces.  Hence,  we  need  to  account  for  this  lack 
of  continuity  when  designing  a  discretization  scheme  for  the  steady  Navier-Stokes 
equations.  We  weakly  enforce  tangential  continuity  between  adjacent  patches  using 
a  combination  of  upwinding  and  the  symmetric  interior  penalty  method  [1,  15,  39]. 

We  now  establish  some  preliminary  notations.  Let  ICh{^i)  and  Jbi(f^i)  denote  the 
sets  of  physical  mesh  elements  and  faces  associated  with  patch  flj.  We  denote  the 
global  set  of  mesh  elements  as  JCh  and  the  global  set  of  mesh  faces  as  Jbi-  As  in  the 
single  patch  setting,  we  define  the  boundary  mesh  to  be 

Th  =  {F  e  Fh{^i),i  =  I, . . .  ,np  :  F  G  (912}  ,  (94) 

and  we  dehne  the  interface  mesh  to  be 

Ih  =  {F  e  =  1, . . .  ,np  :  F  G  7^  j  and  F  ^Th}-  (95) 

For  each  face  F  G  X/j  belonging  to  the  interface  mesh,  there  exist  two  unique  adjacent 
elements  F+,  K~  G  /C^  such  that  F  G  dK^  and  F  G  dK~ .  We  dehne  for  such  a  face 
the  mesh  size 

hp  :=  -  {hK+  +  hx-)  ■  (96) 

Let  (j)  be  an  arbitrary  scalar-,  vector-,  or  matrix-valued  piecewise  smooth  function, 
and  let  us  denote  by  0+  and  0“  the  traces  of  0  on  F  as  taken  from  within  the  interior 
of  K'^  and  K~  respectively.  We  dehne  the  mean  value  of  0  at  x  G  F  as 

}  (G  +  r) .  (97) 
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Further,  for  a  generic  multiplication  operator  O,  we  define  the  jump  of  ^On  at  x  G  -F 
as 

10  0  n]  :=  0+  0  nx+  +  0“  0  n^-  (98) 

where  nj^+/-  denotes  the  outward  facing  normal  on  the  boundary  dK^I~  of  element 
K+/-. 

With  the  above  notation  established,  let  us  dehne  the  following  bilinear  form; 


v)  =  ^  (2z/V"w,  vw; 


(L2(ni))dxd 


i=l 


(§VW}}  :  |w  0  n]  +  {{V"w}  :  |v  0  n])  ds 


F&Xh' 


|w  0  n]  :  |v  0  n|  )  ds 

hjp 


Feih 

[  2z/ (  ((VW)n)  ■  w+ ((V"w)n)  ■  V 
Fer^  ^ 


a 


pen 


hi 


w  ■  V  ]  ds.  (99) 


Above,  Open  >  0  denotes  the  same  positive  penalty  constant  as  before.  We  note  that 
the  above  bilinear  form  is  coercive.  Now,  for  hlj  an  arbitrary  patch,  let  n*  denote  the 
outward-facing  normal  with  respect  to  dVLi.  We  dehne  the  upwind  form 

Up 

4(w,  x;  v)  =  ^  -  (w  0  X,  Vv)(^2(f^^))dxd 
i=l 

+  ^  /  -  (w  ■  rij  |w  ■  n*!)  u  ■  vds 

JdQi\dU  ^ 

f  1 

+  ^  /  -  (w  ■  rij  -  |w  ■  riil)  W  ■  vds 

Jdni\dQ  2 


for  w,  X,  V  G  Ho(div;  hi)  where  u®  is  the  trace  of  u  taken  from  the  exterior  of  hi*.  We 
would  like  to  remark  that  the  above  form  is  nonlinear  in  w.  Furthermore,  the  form 
satishes  the  semi-coercivity  result  c)j(w,v;v)  >  0,  and  c)j(w,x;v)  =  0  for  constant 
vector-valued  functions  v  :  hi  — )■  Hence,  the  upwind  form  is  conservative. 

With  all  the  preceding  terminology  dehned,  our  discrete  multi-patch  formulation 
reads  as  follows. 


f  Find  Uh  G  Vo,/i  and  ph  G  Qo,h  such  that 


{MP){ 


kl{uh,  v/i)  +  c*{uh,  ut;  v/i)  -  b{ph,  v/^)  -h  b{qh,  u^)  =  (f,  v;i)L2(f2) 

(100) 

for  all  \h  e  Vo,h  and  qh  G  Qo,h- 
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As  in  the  single  patch  setting,  the  discrete  formulation  detailed  above  returns  a  point- 
wise  divergence-free  velocity  field.  Furthermore,  the  discrete  formulation  is  consistent. 
However,  we  do  not  yet  have  a  convergence  analysis  available.  We  anticipate  such  a 
convergence  analysis  will  take  new  theoretical  developments.  Nonetheless,  we  have 
utilized  the  above  formulation  in  practice  and  observed  it  produces  optimal  conver¬ 
gence  rates  for  both  velocity  and  pressure  helds. 


8  Numerical  Verification  of  Convergence  Estimates 

In  this  section,  we  numerically  verify  our  convergence  estimates  using  a  collection  of 
problems  with  exact  solutions.  Throughout,  we  choose  Nitsche’s  penalty  constant  as 

Cpen  =  5(/i;  -|-  1) 

which  we  have  found  to  be  sufficiently  large  in  order  to  ensure  numerical  stability. 
Additionally,  we  employ  uniform  parametric  meshes,  linear  parametric  mappings,  and 
B-spline  spaces  of  maximal  continuity. 

8.1  Two-dimensional  Manufactured  Solution 

As  a  hrst  numerical  experiment,  we  consider  a  two-dimensional  manufactured  vortex 
solution  that  was  originally  presented  in  [8].  Let 

Vt  =  (0,1)2 


and 

f  =  V  ■  (u  (g)  u)  -  V  ■  (2z/V®u)  +  Vp 

with 

__  2e^{—l  +  xYx^{y‘^  —  y){—l  +  2y) 

^  (— e^(— 1  -f  x)x{—2  +  a;(3  -|-  a;))(— 1  -|-  y)‘^y‘^) 

and 

p  =  (-424 +  156e+(p2-y)(-456  +  e^(456  +  a;2(228  -  5(p2 -!/))  + 

2a;(-228  +  {y^  -  y))  +  2x\-^Q  +  {y^  -  y))  +  x\l2  +  {y^  -  y))))). 

Homogeneous  boundary  conditions  are  applied  along  the  boundary  dVL,  and  the  condi¬ 
tion  f^pdx  =  0  is  enforced.  A  solution  to  the  steady  Navier-Stokes  problem  with  the 
prescribed  forcing  is  then  clearly  (u,p)  =  (u,p),  and  this  solution  is  unique  provided 
a  smallness  condition  is  satished.  The  streamlines  and  pressure  contours  associated 
with  the  solution  are  plotted  in  Figure  2. 

To  conhrm  our  theoretically  derived  error  estimates,  we  have  computed  conver¬ 
gence  rates  for  divergence-conforming  B-spline  discretizations  of  varying  mesh  size  and 
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X 


(a) 


(b) 

Figure  2:  Vortex  manufactured  solution  in  2-D:  (a)  Flow  velocity  streamlines,  (b) 
Pressure  contours. 
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Table  1:  Steady  vortex  flow  convergence  rates  in  2-D:  Re  =  1 


Polynomial  degree  k'  =  1 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  UhWh 

5.48e-2 

2.80e-2 

1.40e-2 

7.00e-3 

3.50e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

5.48e-2 

2.80e-2 

1.40e-2 

7.00e-3 

3.50e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

1  U  —  U/i| 

2.77e-3 

8.16e-4 

2.28e-4 

6.10e-5 

1.58e-5 

order 

- 

1.76 

1.84 

1.90 

1.95 

\\p-Ph\\LHn) 

5.04e-3 

1.38e-3 

3.49e-4 

8.72e-5 

2.18e-5 

order 

- 

1.87 

1.98 

2.00 

2.00 

Polynomial  degree  k'  =  2 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  U/jII/j 

9.71e-3 

2.33e-3 

5.68e-4 

1.40e-4 

3.48e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

U  —  U/i  Hi(0) 

9.70e-3 

2.33e-3 

5.68e-4 

1.40e-4 

3.48e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

1  u  ~  L2(0) 

2.94e-4 

3.84e-5 

5.03e-6 

6.47e-7 

8.21e-8 

order 

- 

2.94 

2.93 

2.96 

2.98 

\\p-Ph\\LHn) 

1.08e-3 

1.12e-4 

1.17e-5 

1.19e-6 

1.27e-7 

order 

- 

3.40 

3.26 

3.30 

3.23 

Polynomial  degree  k'  =  3 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  U/ill/, 

9.86e-4 

1.28e-4 

1.66e-5 

2.13e-6 

2.72e-7 

order 

- 

2.95 

2.95 

2.96 

2.97 

9.83e-4 

1.28e-4 

1.65e-5 

2.10e-6 

2.66e-7 

order 

- 

2.94 

2.96 

2.97 

2.98 

1  U  —  U/il 

3.05e-5 

2.34e-6 

1.59e-7 

1.03e-8 

6.55e-10 

order 

- 

3.70 

3.88 

3.95 

3.98 

\\p-Ph\\L^n) 

l.lOe-4 

5.64e-6 

3.45e-7 

2.19e-8 

1.39e-9 

order 

- 

4.29 

4.03 

3.98 

3.98 
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Table  2:  Robustness  of  2-D  divergence-free  B-spline  discretizations  for  increasing  Re 


Polynomial  degree  k'  =  1,  h  =  1/1& 


Re 

0 

1 

10 

100 

1000 

10000 

||U  -  Uftll,, 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

u  —  U/j 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

1.40e-2 

1  “  U/ii  l2(0) 

2.28e-4 

2.28e-4 

2.28e-4 

2.28e-4 

2.28e-4 

2.28e-4 

\\p-Ph\\L^(Q) 

3.49e-4 

3.49e-4 

1.98e-4 

1.96e-4 

1.96e-4 

1.96e-4 

Polynomial  degree  k'  =  2,  h  =  1/1& 


Re 

0 

1 

10 

100 

1000 

10000 

1  u 

-  Uftll,, 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

u  — 

U/i  HpO) 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

5.68e-4 

1  u  — 

U/ii  L2(0) 

5.03e-6 

5.03e-6 

5.03e-6 

5.03e-6 

5.03e-6 

5.03e-6 

lb- 

P/ll  L2(Q) 

1.17e-5 

1.17e-5 

6.50e-6 

6.42e-4 

6.42e-6 

6.42e-6 

Polynomial  degree  /c'  =  3,  /i  =  l/16 


Re 

0 

1 

10 

100 

1000 

10000 

||U  -  Uftll,, 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

b  —  Uh  HpO) 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

1.66e-5 

lb  —  u/il  l2(o) 

1.59e-7 

1.59e-7 

1.59e-7 

1.59e-7 

1.59e-7 

1.59e-7 

\\p-Ph\\LHn) 

3.45e-7 

3.45e-7 

3.19e-7 

3.19e-7 

3.19e-7 

3.19e-7 

Table  3:  Instability  of  conservative  2-D  Taylor-Hood  discretizations  for  increasing  Re 


Q2/Q1  velocity /pressure  pair,  h  =  1/16 


Re 

0 

1 

10 

100 

1000 

10000 

b  —  u/i  Hpn) 

6.78e-4 

6.78e-4 

7.11e-4 

2.26e-3 

2.16e-2 

2.16e-l 

lb  —  u/il  L2(n) 

6.54e-6 

6.54e-6 

6.79e-6 

1.97e-5 

1.86e-4 

2.35e-3 

\\p-Ph\\L^(Q) 

1.96e-4 

1.96e-4 

1.96e-4 

1.96e-4 

1.96e-4 

1.96e-4 

35 


polynomial  degree.  Furthermore,  we  have  computed  rates  for  a  variety  of  Reynolds 
numbers 


V 


where  U  is  a  velocity  scale  parameter  and  L  is  a  length  scale  parameter.  We  take  U  = 
1  and  L  =  1.  In  Table  1,  we  have  plotted  our  numerically  computed  convergence  rates 
for  Re  =  1.  First,  note  that  our  theoretically  derived  error  estimates  are  conhrmed. 
Second,  note  that  the  L^-norm  of  the  pressure  error  optimally  converges  like  0{h^ 
which  is  an  improvement  over  our  theoretically  derived  estimate.  Third,  note  that 
the  results  in  Table  1  are  identical  to  the  Stokes  results  appearing  in  Table  2  of  [19]. 
Hence,  the  introduction  of  convection  has  not  affected  our  numerical  error. 

While  our  theoretically  derived  error  estimates  only  cover  flows  with  “small” 
Reynolds  number  (and  indeed  uniqueness  is  only  guaranteed  under  a  smallness  condi¬ 
tion),  we  have  investigated  the  effectiveness  of  our  discretization  for  larger  Reynolds 
number  flows  and  the  two-dimensional  forcing  prescribed  above.  To  compute  these 
flow  solutions,  we  used  a  Newton-Raphson  nonlinear  solver  in  conjunction  with  con¬ 
tinuation.  The  results  of  this  investigation  are  included  in  Table  2  for  meshes  with 
16  X  16  elements.  Note  that  the  velocity  numerical  errors  in  the  tables  appear  in¬ 
dependent  of  the  Reynolds  number.  Moreover,  the  pressure  numerical  errors  seem 
to  improve  with  increasing  Reynolds  number.  Hence,  we  are  able  to  recover  the  de¬ 
sired  flow  solution  (u,  p)  for  large  Reynolds  numbers.  This  observation  attests  to  the 
enhanced  stability  properties  of  our  discretization  for  nonlinear  flow  problems  even 
in  the  absence  of  any  external  stabilization  mechanisms.  To  contrast  our  methodol¬ 
ogy  with  standard  mixed  hnite  element  discretizations,  we  have  repeated  the  above 
computations  for  conservative  Taylor-Hood  [24]  hnite  element  approximations,  us¬ 
ing  again  a  continuation  method  in  conjunction  with  a  Newton-Raphson  nonlinear 
solver.  The  results  of  these  computations  are  included  in  Table  3.  Note  that  while  the 
pressure  error  is  well-behaved,  the  velocity  error  diverges  with  increasing  Reynolds 
number.  We  believe  that  this  divergence  is  inherently  tied  to  the  fact  the  Taylor-Hood 
approximations  only  satisfy  the  divergence-free  constraint  approximately. 

8.2  Three-dimensional  Manufactured  Solution 

As  a  second  numerical  experiment,  we  consider  a  three-dimensional  manufactured 
solution  representing  a  single  vortical  hlament.  Let 

Q  =  (0,1)3 


and 

f  =  V  ■  (u  (g)  u)  -  V  ■  (2z/V*u)  -h  Vp 


with 


u  =  curl0. 


36 


Figure  3:  Vortex  manufactured  solution  in  3-D:  Flow  velocity  streamlines  colored  by 
velocity  magnitude. 


x{x  —  l)y‘^{ri  —  —  1)^ 

0 

x^{x  —  —  V)^z{z  —  1) 


and 


Again,  homogeneous  boundary  conditions  are  applied  along  the  boundary  dfl,  and 
the  pressure  is  enforced  to  satisfy  J^pd^s.  =  0.  A  solution  to  the  steady  Navier-Stokes 
equations  with  the  prescribed  forcing  is  then  clearly  (u,  p)  =  (u,  p),  and  for  sufficiently 
small  data,  this  solution  is  unique.  Streamlines  associated  with  the  exact  solution  are 
plotted  in  Figure  3. 

For  Re  =  1  flow  (with  Re  =  1/z^),  we  have  computed  convergence  rates  for 
divergence-conforming  B-spline  discretizations  of  varying  mesh  size  and  polynomial 
degree.  In  Table  4,  we  have  listed  our  numerically  computed  convergence  rates. 
First,  note  that  our  theoretically  derived  error  estimates  are  confirmed.  Second,  note 
that  the  L^-norm  of  the  pressure  error  optimally  converges  like  which  is 

an  improvement  over  our  theoretically  derived  estimate.  Third,  note  that,  as  in  the 
two-dimensional  setting,  the  results  in  Table  4  are  identical  to  the  Stokes  results 
appearing  in  Table  5  of  [19].  Hence,  the  introduction  of  convection  has  not  affected 
our  numerical  error. 

Again,  while  our  theoretically  derived  error  estimates  only  cover  flows  with  “small” 
Reynolds  number,  we  have  investigated  the  effectiveness  of  our  discretization  for  larger 
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Table  4:  Steady  vortex  flow  convergence  rates  in  3-D;  Re  =  1. 


Polynomial  degree  k'  =  1 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh\\h 

1.83e-2 

8.98e-3 

4.18e-3 

1.99e-3 

9.62e-4 

order 

- 

1.03 

1.10 

1.07 

1.05 

1.51e-2 

7.64e-3 

3.77e-3 

1.87e-3 

9.34e-4 

order 

- 

0.98 

1.02 

1.01 

1.00 

1.35e-3 

3.68-4 

1.03e-4 

2.81e-5 

7.40e-6 

order 

- 

1.88 

1.84 

1.87 

1.93 

\\p-Ph\\LHn) 

5.41e-2 

1.48e-2 

3.58e-3 

8.85e-4 

2.26e-4 

order 

- 

1.87 

2.05 

2.02 

1.97 

Polynomial  degree  k'  =  2 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

||U  -  Uh\\h 

6.50e-3 

1.54e-3 

4.10e-4 

9.51e-5 

2.15e-5 

order 

- 

2.08 

1.91 

2.11 

2.15 

U  —  U/i 

3.71e-3 

9.90e-4 

2.79e-4 

6.59e-5 

1.50e-5 

order 

- 

1.91 

1.83 

2.08 

2.14 

D(o) 

1.97e-4 

4.25e-5 

7.38e-6 

8.67e-7 

9.18e-8 

order 

- 

2.21 

2.53 

3.09 

3.23 

\\p-Ph\\LHn) 

1.50e-2 

1.59e-3 

2.00e-4 

2.56e-5 

3.26e-6 

order 

- 

3.24 

2.99 

2.97 

2.97 

Table  5:  Robustness  of  3-D  divergence-free  B-spline  discretizations  for  increasing  Re 


Polynomial  degree  k'  =  1,  h  =  1/16 


Re 

0 

1 

10 

100 

1000 

10000 

i|u  -  Uftll/j 

1.99e-3 

1.99e-3 

1.99e-3 

1.99e-3 

1.99e-3 

1.99e-3 

u  —  U/i  HpO) 

1.87e-3 

1.87e-3 

1.87e-3 

1.87e-3 

1.87e-3 

1.87e-3 

2.81e-5 

2.81e-5 

2.81e-5 

2.81e-5 

2.81e-5 

2.81e-5 

\\P-Ph\\mn) 

8.84e-4 

8.84e-4 

8.84e-4 

8.84e-4 

8.84e-4 

8.84e-4 

Polynomial  degree  k'  =  2,  h  =  1/16 


Re 

0 

1 

10 

100 

1000 

10000 

1  u 

-  Uftll/i 

9.51e-5 

9.51e-5 

9.51e-5 

9.51e-5 

9.51e-5 

9.51e-5 

u  — 

U/i  HpO) 

6.59e-5 

6.59e-5 

6.59e-5 

6.59e-5 

6.59e-5 

6.59e-5 

1  u  — 

L2(0) 

8.67e-7 

8.67e-7 

8.67e-7 

8.67e-7 

8.67e-7 

8.67e-7 

lb- 

P/i|  L2(0) 

2.56e-5 

2.56e-5 

2.56e-5 

2.56e-5 

2.56e-5 

2.56e-5 
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Reynolds  number  flows  and  the  three-dimensional  forcing  prescribed  above.  We  used 
a  Newton-Raphson  nonlinear  solver  in  conjunction  with  continuation.  The  results 
of  this  investigation  are  included  in  Table  5  for  meshes  with  16  x  16  x  16  elements. 
Note  that  the  velocity  and  pressure  numerical  errors  appearing  in  the  tables  are 
independent  of  the  Reynolds  number.  Hence,  we  are  able  to  recover  the  desired  flow 
solution  (u,p)  for  large  Reynolds  numbers. 

8.3  Kovasznay  Flow 

As  a  third  numerical  experiment,  we  consider  Kovasznay  flow.  Kovasznay  flow  refers 
to  the  flow  behind  an  infinite  two-dimensional  grid,  and  it  is  often  utilized  as  a 
convergence  test  for  Navier-Stokes  discretizations.  The  flow  pattern  is  periodic  and 
can  be  analytically  derived  [27].  Indeed,  denoting  Re  =  the  flow  solution  satisfies 

_  1  —  cos(27r|/) 

^  .  ^e^"'cos(27r|/). 


and 


where 

Re  /Re^  ~ 

A-  — -y— +47r. 

The  streamlines  and  pressure  contours  associated  with  the  Kovasznay  flow  solution 
at  Re  =  40  are  plotted  in  Figure  4.  Note  that  the  streamlines  closely  resemble  the 
streamlines  associated  with  steady  flow  behind  a  cylinder. 

To  solve  the  Kovasznay  flow  problem  numerically,  we  restrict  ourselves  to  the 
domain  H  =  (0, 1)  x  (—1/2, 1/2).  On  the  left,  bottom,  and  top  sides  of  the  domain, 
we  enforce  Dirichlet  boundary  conditions.  On  the  right  side  of  the  domain,  we  enforce 
a  traction  boundary  condition.  On  the  interior  of  the  domain,  we  apply  the  usual 
forcing 

f  =  V  ■  (u  0  u)  -  V  ■  (2z/V*u)  -h  Vp. 

We  have  computed  Re  =  40  flow  convergence  rates  for  a  variety  of  divergence- 
conforming  B-spline  discretizations.  These  rates  are  summarized  in  Table  6.  Note 
that  the  convergence  rates  are  approaching  the  optimal  rates  as  h  — )■  0  for  both 
the  velocity  and  pressure  held.  However,  it  is  apparent  that  even  for  h  =  1/64  our 
computations  still  lie  in  the  pre-asymptotic  range. 


8.4  Cylindrical  Couette  Flow 

As  a  final  convergence-rate  experiment,  we  consider  cylindrical  Couette  flow.  Couette 
flow  is  often  used  as  a  “sanity  check”  for  Navier-Stokes  discretizations.  Cylindrical 
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(a) 


0.5| 


0.2 


0.1 


>.  0 


-0.1 
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-0.5  ^ 
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0.35 


0.25 


0.2 


0.15 


(b) 


Figure  4:  Steady  Kovasznay  flow:  (a)  Streamlines  for  Re  =  40,  (b)  Pressure  contours 
for  i?e  =  40. 


40 


Table  6:  Steady  Kovasznay  flow  convergence  rates:  Re  =  40 


Polynomial  degree  k'  =  1 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

u  —  U/i  h1(Q) 

1.39e0 

7.31e-l 

3.69e-l 

1.84e-l 

9.19e-2 

order 

- 

0.93 

0.99 

1.00 

1.00 

5.31e-2 

1.98e-2 

6.78e-3 

2.15e-3 

6.34e-4 

order 

- 

1.43 

1.54 

1.66 

1.76 

\\p-Ph\\LHn) 

3.98e-2 

1.49e-2 

4.73e-3 

1.35e-3 

3.75e-4 

order 

- 

1.42 

1.65 

1.81 

1.85 

Polynomial  degree  k'  =  2 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

U  —  U/i  HpO) 

4.59e-l 

1.17e-l 

2.78e-2 

6.69e-3 

1.64e-3 

order 

- 

1.97 

2.08 

2.05 

2.03 

1.44e-2 

1.96e-3 

2.41e-4 

3.04e-5 

2.83e-6 

order 

- 

2.88 

3.02 

2.99 

2.99 

\\p-Ph\\LHn) 

1.65e-2 

3.55e-3 

5.14e-4 

7.05e-5 

9.56e-6 

order 

- 

2.22 

2.79 

2.87 

2.88 

Polynomial  degree  k'  =  3 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

U  —  U/i 

1.29e-l 

1.52e-2 

1.94e-3 

2.55e-4 

3.31e-5 

order 

- 

3.08 

2.98 

2.92 

2.95 

l2(0) 

2.95e-3 

1.97e-4 

1.64e-5 

1.20e-6 

8.51e-8 

order 

- 

3.91 

3.59 

3.77 

3.81 

Wp-PhWL-^in) 

5.59e-3 

6.48e-4 

5.75e-5 

5.28e-6 

4.00e-7 

order 

- 

3.11 

3.49 

3.45 

3.72 
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Couette  flow  is  a  more  realistic  problem  which  describes  the  flow  between  two  concen¬ 
tric  rotating  cylinders.  Here,  we  consider  flow  between  a  fixed  outer  cylinder  and  a 
rotating  inner  cylinder.  The  problem  setup  is  illustrated  in  Figure  5(a).  No  external 
forcing  is  applied.  The  velocity  field  for  this  flow  assumes  the  form 


u  = 


where 


U0{r)  sin(6*) 
Ud{r)  cos{9) 

ue{r)  =  Ar 


(r,  9)  are  polar  coordinates  with  respect  to  the  center  of  the  cylinders,  and 


A  =  -Q, 


<52 


52’ 


D  _  O  O  _  ^ 

ID  —  /-,  ro\  ?  —  1  ^  — 

(1  -  C>2)  Tin  Tout 


We  have  depicted  this  velocity  held  in  Figure  5(b).  The  pressure  held  for  cylindrical 
Couette  how  is  axisymmetric  and  satishes  the  diherential  equation 

(101) 


dp{r) 


ue[r) 


dr  r 

The  above  diherential  equation  coupled  with  the  constraint 


0 


determines  the  pressure  held  uniquely.  In  what  follows,  we  assume  Tin  =  1,  Tout  =  2, 
and  U  =  1. 

We  have  computed  convergence  rates  for  a  variety  of  divergence-conforming  B- 
spline  discretizations  and  for  Re  =  40.  To  represent  the  annular  domain  in  our 
computations,  we  employed  the  polar  mapping 


F(6,6) 


{{Tout  -  rin)i2  +  Tin)  sin(27r^i) 
{{rout  -  nn)6  +  rin)  cos(27r^i) 


,V(ei,6)e  (0,1)2 


(102) 


and  periodic  B-splines  of  maximal  continuity  in  the  .^i-direction  (see  Section  2  of  [17]). 
It  should  be  emphasized  that  we  do  not  use  the  polar  form  of  the  steady  Navier-Stokes 
equations.  Rather,  we  utilize  the  polar  mapping  to  dehne  our  divergence-conforming 
B-splines  in  physical  space  and  then  employ  the  Cartesian-based  variational  formu¬ 
lation  discussed  in  this  chapter.  The  results  of  our  computations  are  summarized 
in  Table  7.  Note  from  the  table  that  all  of  our  theoretically  derived  error  estimates 
are  conhrmed.  Furthermore,  we  have  obtained  axisymmetric  velocity  helds  with  zero 
radial  components,  and  the  discrete  pressure  held  converges  at  optimal  order.  We 
have  additionally  simulated  the  cylindrical  Couette  how  problem  using  a  multi-patch 
NURBS  construction  (see  Subsection  8.4  of  [19])  and  obtained  identical  discrete  ve¬ 
locity  helds  and  slightly  dihering  (yet  still  optimally  convergent)  discrete  pressure 
helds. 
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(b) 


Figure  5:  Cylindrical  Couette  flow:  (a)  Problem  setup,  (b)  Flow  velocity  arrows. 
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h/ho  =  1/8 

Figure  6:  Cylindrical  Couette  flow:  Sequence  of  polar  meshes. 
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Table  7:  Cylindrical  Couette  flow  convergence  rates:  Re  =  40 


Polynomial  degree  k'  =  1 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

||U  -  Uh\\h 

5.42e-l 

2.63e-l 

1.26e-l 

6.12e-2 

2.99e-2 

order 

- 

1.04 

1.06 

1.04 

1.03 

U  —  U/i 

4.48e-l 

2.32e-l 

1.17e-l 

5.86e-2 

2.93e-2 

order 

- 

0.95 

0.99 

1.00 

1.00 

1^1  “  l2(0) 

5.00e-2 

1.53-2 
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Figure  7:  Lid-driven  flow  in  a  two-dimensional  cavity:  Problem  setup. 

9  Benchmark  Problems 

In  this  section,  we  investigate  the  effectiveness  of  our  methodology  as  applied  to 
two  benchmark  problems:  two-dimensional  lid-driven  cavity  flow  and  confined  jet 
impingement.  As  previously,  we  choose  Nitsche’s  penalty  constant  as  Open  =  5(A;'-|-1) 
for  all  of  the  following  examples,  and  we  employ  uniform  parametric  meshes,  linear 
parametric  mappings,  and  B-spline  spaces  of  maximal  continuity. 

9.1  Two-dimensional  Lid-driven  Cavity  Flow 

Two-dimensional  lid-driven  cavity  flow  is  widely  considered  to  be  one  of  the  classical 
test  problems  for  the  validation  of  numerical  discretizations  for  incompressible  flow 
simulation.  In  the  presence  of  increasing  Reynolds  number,  lid-driven  cavity  flow 
loses  its  symmetry  and  eventually  becomes  unsteady.  The  problem  setup  for  lid- 
driven  cavity  flow  is  shown  in  Figure  7.  For  the  computations  here,  H  and  U  are 
defined  to  be  1.  The  Reynolds  number  associated  with  the  flow  is  defined  to  be 

V 

The  left,  right,  and  bottom  sides  of  the  cavity  are  fixed  no-slip  walls  while  the  top 
side  of  the  cavity  is  a  wall  which  slides  to  the  right  with  velocity  magnitude  U.  The 
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Figure  8:  Lid-driven  cavity  flow:  Streamlines  for  Re  =  100. 

forcing  f  is  defined  to  be  zero.  The  pressure  and  stress  fields  associated  with  this  flow 
experience  corner  singularities  which  impede  the  convergence  of  numerical  methods 
and  expose  unstable  velocity/pressure  pairs.  Furthermore,  the  velocity  boundary 
condition  is  discontinuous  at  the  upper  left  and  right  corners  of  the  domain.  In 
classical  finite  element  analysis,  this  boundary  condition  must  be  replaced  with  an 
approximate  “smooth”  boundary  condition.  Here,  as  the  prescribed  slip  boundary 
condition  is  weakly  enforced,  we  do  not  have  to  invoke  such  an  approximation  and 
can  instead  directly  utilize  the  discontinuous  boundary  condition. 

We  have  numerically  simulated  lid-driven  cavity  flow  at  Re  =  100,  Re  =  400,  and 
Re  =  1000  using  a  variety  of  divergence-conforming  B-spline  discretizations.  Stream¬ 
lines  computed  using  a  mesh  of  degree  k'  =  1  and  size  h  =  1/128  are  plotted  in  Figures 
8,  9,  and  10  for  Re  =  100,  Re  =  400,  and  Re  =  1000,  respectively.  These  computed 
streamlines  are  visually  indistinguishable  from  benchmark  streamlines  appearing  in 
the  literature. 

We  have  compared  our  numerical  results  with  the  classical  benchmark  results  of 
Ghia  et  al.  [21]  for  three  different  polynomial  degrees  {k'  =  1,  k'  =  2,  k'  =  3), 
three  different  mesh  sizes  {h  =  1/32,  h  =  1/64,  and  h  =  1/128),  and  the  selected 
Reynolds  numbers.  The  results  of  Ghia  were  obtained  using  a  second-order  upwind 
finite  difference  method  on  a  stretched  mesh  with  129^  grid  points.  In  Figures  11(a), 
12(a),  and  13(a),  we  have  compared  our  k'  =  1,  Re  =  100  values  for  the  horizontal 
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Figure  9:  Lid-driven  cavity  flow:  Streamlines  for  Re  =  400. 

component  of  the  velocity  field  along  the  vertical  center  line  with  those  of  Ghia, 
and  in  Figures  11(b),  12(b),  and  13(b),  we  have  compared  our  k'  =  1,  Re  =  100 
values  for  the  vertical  component  of  the  velocity  held  along  the  horizontal  center 
line.  Note  that  our  centerline  velocity  results  are  roughly  independent  of  the  mesh 
size.  Second,  note  that  while  our  results  ostensibly  match  those  of  Ghia,  there  are 
discernible  quantitative  differences  between  the  two  sets  of  results  even  for  our  finest 
mesh.  To  investigate  these  differences  further,  we  have  employed  a  set  of  converged 
pseudospectral  results  that  were  obtained  via  a  subtraction  method  to  remove  the 
leading  terms  of  the  corner  singularities  [5].  In  Table  8,  we  have  compared  our 
centerline  results  with  these  converged  results  as  well  as  the  results  of  Ghia.  Note 
that  our  results  are  much  more  accurate  than  the  results  of  Ghia  for  all  mesh  sizes 
and  polynomial  degrees.  Furthermore,  our  results  for  k'  =  2  and  /c'  =  3  on  a  64  x  64 
element  mesh  are  indistinguishable  from  the  converged  pseudospectral  results.  This 
attests  to  the  effectiveness  of  our  discretization  with  increasing  k' ,  even  in  the  presence 
of  singularities.  We  should  also  again  remark  that  these  results  were  obtained  without 
stretched  meshes  and  without  stabilization.  We  believe  that  the  combination  of  exact 
mass  conservation  and  weak  enforcement  of  the  no-slip  condition  plays  a  pivotal  role 
in  the  enhanced  accuracy  of  our  discretization  scheme. 

In  Figures  14(a),  15(a),  and  16(a),  we  have  compared  our  k'  =  1,  Re  =  400  values 
for  the  horizontal  component  of  the  velocity  field  along  the  vertical  center  line  with 
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Figure  10:  Lid-driven  cavity  flow:  Streamlines  for  Re  =  1000. 

those  of  Ghia,  and  in  Figures  14(b),  15(b),  and  16(b),  we  have  compared  our  k'  =  1, 
Re  =  400  values  for  the  vertical  component  of  the  velocity  held  along  the  horizontal 
center  line.  Note  that  our  Re  =  400  results  match  the  results  of  Ghia  better  than 
the  Re  =  100  results.  In  fact,  our  Re  =  400  results  are  nearly  indistinguishable  from 
the  benchmark  results  on  the  finest  mesh.  This  is  encouraging  as  the  Ghia  Re  =  400 
results  match  well  with  highly-accurate  extrapolated  results  in  the  literature  [34]. 
Moreover,  our  results  for  the  coarsest  mesh  match  the  results  of  Ghia  better  than  any 
comparable  second-order  results  we  have  seen  in  the  literature.  In  Table  8,  we  have 
compared  our  centerline  results  with  the  results  of  Ghia  along  with  our  results  for 
k'  =  2  and  /c'  =  3  on  a  64  x  64  element  mesh.  Note  that  our  results  appear  to  converge 
quickly  with  increasing  k',  despite  the  increased  smoothness  of  our  discrete  spaces. 
By  using  our  k'  =  3  results  as  a  benchmark,  we  see  our  k'  =  1  results  are  considerably 
more  accurate  than  Ghia’s  results  for  h  <  1/64.  Unfortunately,  no  pseudospectral 
results  are  available  to  use  as  comparison. 

In  Figures  17(a),  18(a),  and  19(a),  we  have  compared  our  k'  =  1,  Re  =  1000  values 
for  the  horizontal  component  of  the  velocity  held  along  the  vertical  center  line  with 
those  of  Ghia,  and  in  Figures  17(b),  18(b),  and  19(b),  we  have  compared  our  k'  =  1, 
Re  =  1000  values  for  the  vertical  component  of  the  velocity  held  along  the  horizontal 
center  line.  This  is  a  more  challenging  test  case  than  either  the  Re  =  100  or  Re  =  400 
flows,  and  the  coarse  32  x  32  element  mesh  is  not  nearly  fine  enough  to  resolve  the 
flow  features.  Nonetheless,  our  results  for  this  coarse  mesh  better  match  the  results 
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Table  8:  Two-dimensional  lid-driven  cavity  flow:  Extrema  of  the  velocity  through  the 
centerlines  of  the  cavity. 


Re  =  100 


Method 

^min 

^max 

^min 

B-spline,  k'  =  1  and  h  =  1/32 

-0.21551 

0.18054 

-0.25472 

B-spline,  k'  =  1  and  h  =  1/64 

-0.21443 

0.17991 

-0.25409 

B-spline,  k'  =  1  and  h  =  1/128 

-0.21414 

0.17966 

-0.25387 

B-spline,  k'  =  2  and  h  =  1/64 

-0.21404 

0.17957 

-0.25379 

B-spline,  k'  =  3  and  h  =  1/64 

-0.21404 

0.17957 

-0.25380 

Pseudospectral  (Ref.  [5]) 

-0.21404 

0.17957 

-0.25380 

Ghia  et  al.  (Ref.  [21]) 

-0.21090 

0.17527 

-0.24533 

Re  = 

=  400 

Method 

^min 

'^raax 

"^min 

B-spline,  k'  =  1  and  h  =  1/32 

-0.33651 

0.31039 

-0.45768 

B-spline,  k'  =  1  and  h  =  1/64 

-0.33150 

0.30605 

-0.45659 

B-spline,  k'  =  1  and  h  =  1/128 

-0.32989 

0.30471 

-0.45470 

B-spline,  k'  =  2  and  h  =  1/64 

-0.32927 

0.30415 

-0.45406 

B-spline,  k'  =  3  and  h  =  1/64 

-0.32925 

0.30413 

-0.45401 

Ghia  et  al.  (Ref.  [21]) 

-0.32376 

0.30203 

-0.44993 

Re  = 

1000 

Method 

'^min 

'^max 

'^min 

B-spline,  k'  =  1  and  h  =  1/32 

-0.40140 

0.39132 

-0.54261 

B-spline,  k'  =  1  and  h  =  1/64 

-0.39399 

0.38229 

-0.53353 

B-spline,  k'  =  1  and  h  =  1/128 

-0.39021 

0.37856 

-0.52884 

B-spline,  k'  =  2  and  h  =  1/64 

-0.38874 

0.37715 

-0.52726 

B-spline,  k'  =  3  and  h  =  1/64 

-0.38857 

0.37698 

-0.52696 

Pseudospectral  (Ref.  [5]) 

-0.38857 

0.37694 

-0.52707 

Ghia  et  al.  (Ref.  [21]) 

-0.38289 

0.37095 

-0.51550 
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Figure  11:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  1,  h  =  1/32,  and  Re  =  100. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  12:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  1,  h  =  1/64,  and  Re  =  100. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  13:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  I,  h  =  1/128,  and  Re  =  100. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  14:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  1,  h  =  1/32,  and  Re  =  400. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  15:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  1,  h  =  1/64,  and  Re  =  400. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  16:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  I,  h  =  1/128,  and  Re  =  400. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  17:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  I,  h  =  1/32,  and  Re  =  1000. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  18:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  I,  h  =  1/64,  and  Re  =  1000. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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Figure  19:  Lid-driven  cavity  flow:  Velocity  fleld  for  k'  =  h  =  1/128,  and  Re  =  1000. 
(a)  Value  of  the  horizontal  component  of  the  velocity  fleld  along  the  vertical  center 
line,  (b)  Value  of  the  vertical  component  of  the  velocity  fleld  along  the  horizontal 
center  line. 
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of  Ghia  than  any  comparable  second-order  results  we  have  seen  in  the  literature.  Our 
results  for  the  64  x  64  and  128  x  128  element  meshes  match  Ghia’s  results  even  better. 
In  Table  8,  we  have  compared  our  centerline  results  with  the  results  of  Ghia  along 
with  a  set  of  converged  pseudospectral  results  with  singularity  subtraction  [5].  As 
expected,  the  results  of  Ghia  are  more  accurate  than  our  coarse  mesh  results,  but  our 
results  for  the  64  x  64  element  mesh  are  comparable  in  accuracy  and  our  results  for 
the  128  X  128  element  mesh  are  considerably  more  accurate.  We  also  computed  results 
for  the  64  x  64  element  mesh  for  polynomial  degrees  k'  =  2  and  k'  =  3  and  found  our 
results  quickly  improved  with  increasing  polynomial  degree.  Indeed,  our  results  for 
k'  =  3  are  nearly  indistinguishable  from  the  converged  pseudospectral  results. 

To  compute  all  of  these  flow  examples,  we  utilized  the  Newton-Raphson  method 
in  conjunction  with  continuation.  Specihcally,  for  the  Re  =  100  simulations,  we 
employed  Newton-Raphson  using  the  results  of  a  corresponding  Stokes  simulation  as 
an  initial  guess.  For  the  Re  =  400  simulations,  we  employed  Newton-Raphson  using 
the  results  of  the  Re  =  100  simulations  as  an  initial  guess,  and  so  on.  We  found 
a  maximum  of  four  Newton-Raphson  steps  were  required  to  achieve  a  sufficiently 
accurate  solution  for  each  nonlinear  solve.  Using  this  procedure,  we  have  been  able 
to  successfully  simulate  flows  upwards  of  Re  =  3200  on  relatively  coarse  meshes.  We 
have  also  been  able  to  reproduce  all  of  the  results  presented  here  by  a  more  expensive 
dynamic  approach  in  which  the  solution  is  evolved  from  an  initial  condition  by  the 
unsteady  Navier-Stokes  equations. 

9.2  Confined  Jet  Impingement 

Impinging  jets  are  commonly  used  in  engineering  applications  due  to  their  enhanced 
heat  and  mass  transfer  characteristics.  Impinging  jets  even  occur  in  blood  vessels  and 
are  believed  to  play  a  role  in  atherogenesis  [20].  Many  discretization  schemes  yield 
a  velocity  held  with  spurious  nonzero  divergence  when  applied  to  the  simulation  of 
incompressible  jet  impingement,  even  at  small  and  moderate  Reynolds  numbers.  As 
our  proposed  discretization  satishes  the  divergence-free  constraint  exactly,  it  does  not 
suffer  from  this  non-physical  behavior. 

The  physical  problem  of  conhned  jet  impingement  is  illustrated  in  Figure  20.  Fluid 
hows  from  one  pipe  into  a  second  pipe  lying  perpendicular  to  the  hrst.  To  simulate 
this  how,  we  use  a  two-dimensional  mathematical  idealization  which  is  illustrated  in 
Figure  21.  Along  the  left  hand  side  of  the  domain,  a  symmetry  condition  is  applied 
as  we  only  model  half  of  the  jet  system.  Along  the  bottom  side  of  the  domain,  no- 
slip  and  no-penetration  boundary  conditions  are  enforced.  Along  the  top  side  of  the 
domain,  two  diherent  sets  of  boundary  conditions  are  applied.  Along  the  hrst  D /2 
length  of  the  top  side,  an  inhow  boundary  condition  is  applied.  Along  the  remainder 
of  the  top  side,  no-slip  and  no-penetration  boundary  conditions  are  enforced.  Along 
the  right  hand  side  of  the  domain,  a  zero-traction  boundary  condition  is  enforced. 
The  height  of  the  domain  is  set  as  H,  and  the  length  of  the  domain,  L,  is  chosen  long 
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Figure  20:  Confined  jet  impingement:  Physical  setup. 
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Figure  21:  Confined  jet  impingement:  Model  problem  description. 
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enough  such  that  the  flow  exiting  the  pipe  is  parallel  to  the  outflow  plane.  For  the 
computations  here,  D  =  H  =  1  and  L  =  8.  The  inflow  speed  is  set  to  be  17  =  1.  The 
Reynolds  number  for  the  flow  is  defined  as  Re  =  U H / v. 

We  have  simulated  Re  =  50  confined  jet  impingement  for  two  divergence-conforming 
B-spline  discretizations  of  polynomial  degree  k'  =  1  corresponding  to  mesh  lengths  of 
h  =  1/16  and  h  =  1/32.  The  resulting  velocity  fields  are  plotted  in  Figures  22  and  23. 
First,  note  that  the  computed  streamlines  for  the  two  meshes  are  indistinguishable. 
Second,  note  that  the  velocity  fields  are  smooth  and  free  from  spurious  oscillations. 
To  compare  our  simulation  results  with  those  of  a  classical  mixed  method,  we  have 
simulated  the  impingement  problem  for  a  conservative  Taylor-Hood  finite  element 
discretization  with  mesh  length  h  =  1/8.  This  discretization  contains  approximately 
the  same  number  of  velocity  degrees  of  freedom  as  the  preceding  coarse  B-spline 
discretization.  The  velocity  held  resulting  from  the  Taylor-Hood  discretization  is 
plotted  in  Figure  24.  Note  that  the  contours  of  the  velocity  held  roughly  match  the 
velocity  contours  resulting  from  the  divergence-conforming  B-spline  discretizations. 
However,  the  velocity  held  resulting  from  the  Taylor-Hood  discretization  exhibits  sig¬ 
nificant  non-physical  compressibility.  We  have  conducted  numerical  investigations  of 
confined  jet  impingement  at  other  Reynolds  numbers  and  obtained  similar  results  to 
those  reported  here. 

10  Conclusions 

In  this  paper,  divergence-conforming  B-spline  discretizations  have  been  presented 
for  the  steady  Navier-Stokes  equations  utilizing  a  variational  formulation  written  in 
conservation  form.  Tangential  velocity  boundary  conditions  are  implemented  weakly 
using  Nitsche’s  method.  The  formulation  yields  velocity  fields  that  are  pointwise 
divergence-free  on  general  mapped  geometries.  A  collection  of  stability  and  error 
estimates  have  been  derived  for  flows  subject  to  a  smallness  condition,  and  these 
theoretical  results  have  been  confirmed  and  supplemented  by  numerical  simulations 
of  problems  with  known  analytical  solutions.  In  fact,  these  numerical  simulations 
have  revealed  that  the  proposed  discretizations  are  robust  with  respect  to  Reynolds 
number,  a  property  not  shared  by  classical  finite  elements  such  as  the  Taylor-Hood 
element.  The  new  discretizations  have  also  been  applied  to  a  selection  of  benchmark 
problems  where  the  advantages  of  the  new  methodology  over  classical  methods  have 
been  highlighted.  In  future  work,  we  plan  to  extend  the  present  developments  to  the 
unsteady  Navier-Stokes  equations  and  generalize  to  locally-refined  meshes. 
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